Additional file 2.

Title: Complete analysis scripts.

Description: This file contains the complete computational pipeline used in this study,
concatenated in execution order. Every script is preceded by a delimiter line of the form
"# FILE: <script name>". The pipeline covers (i) construction of the protein-coding
glycosylation-related gene (GRG) list from the Gene Ontology and UniProt branches,
(ii) microarray pre-processing and quantile normalisation, (iii) differential expression
analysis and functional enrichment, (iv) feature selection by univariate logistic
filtering, random forest ranking and LASSO penalisation, (v) model evaluation with
leak-free repeated cross-validation, decision curve analysis and calibration,
(vi) external validation in GSE30528 and GSE30529, (vii) immune deconvolution using the
CIBERSORT nu-support-vector-regression algorithm with an ssGSEA comparator,
(viii) robustness analyses (bootstrap selection stability, global-factor audit and
scale sensitivity), and (ix) figure and additional-file generation. R and Python package
versions are given in the Methods section of the main manuscript.

Note on paths: machine-specific absolute paths have been generalised. The placeholder
PROJECT_ROOT denotes the root directory of the project on the machine where the analysis
was run; replace it with the local path of the unpacked project before executing the
scripts. Library paths pointing to a personal R library have been written as
"~/software/R_libs". All input data are the public GEO datasets GSE96804, GSE30528 and
GSE30529.

The CIBERSORT R source is redistributed here in its patched form (quantile normalisation
re-implemented in pure R); the original license text is retained in its header.

================================================================================
ANALYSIS SCRIPTS
================================================================================


############################################################
# FILE: CIBERSORT_patched.R
############################################################
# CIBERSORT R script v1.04 (last updated 10-24-2016)

# --- patched quantile normalization (pure R, replaces preprocessCore pthread version) ---
qn_matrix <- function(x) {
  rks <- apply(x, 2, rank, ties.method = "average")
  sm <- rowMeans(apply(x, 2, sort))
  apply(rks, 2, function(r) sm[r])
}
# Note: Signature matrix construction is not currently available;
#   use java version for full functionality.
# Author: Aaron M. Newman, Stanford University (amnewman@stanford.edu)
# Requirements:
#   R v3.0 or later.
# License: https://cibersortx.stanford.edu/CIBERSORT_License.txt

#' Core Algorithm for CIBERSORT Deconvolution
#'
#' @description
#' Performs nu-regression using support vector machines (SVM) to estimate
#' cell type proportions. This is the core computational engine of CIBERSORT,
#' using nu-SVR with linear kernel to decompose mixed gene expression signals.
#'
#' @param X Matrix or data frame containing signature matrix (predictor variables).
#'   Rows are genes, columns are cell types.
#' @param y Numeric vector containing the mixture sample expression (response variable).
#' @param absolute Logical indicating whether to use absolute space for weights.
#'   Default is FALSE (relative proportions).
#' @param abs_method String specifying the method for absolute space weights:
#'   `"sig.score"` or `"no.sumto1"`.
#'
#' @return List containing:
#' \describe{
#'   \item{w}{Estimated cell type weights/proportions}
#'   \item{mix_rmse}{Root mean squared error of the deconvolution}
#'   \item{mix_r}{Correlation coefficient between observed and predicted mixture}
#' }
#'
#' @export
#'
#' @examples
#' # Simulate data
#' set.seed(123)
#' X <- matrix(rnorm(100), nrow = 10)
#' rownames(X) <- paste0("Gene", 1:10)
#' colnames(X) <- paste0("Cell", 1:10)
#' y <- rnorm(10)
#' names(y) <- paste0("Gene", 1:10)
#'
#' # Run core algorithm
#' result <- CoreAlg(X, y, absolute = FALSE, abs_method = "sig.score")
#' if (!is.null(result)) str(result)
CoreAlg <- function(X, y, absolute, abs_method) {
  if (is.null(X)) return(NULL)
  rlang::check_installed("e1071")

  svn_itor <- 3
  nus <- c(0.25, 0.5, 0.75)

  res <- function(i) {
    model <- e1071::svm(X, y,
      type = "nu-regression", kernel = "linear",
      nu = nus[i], scale = FALSE
    )
    model
  }

  if (Sys.info()["sysname"] == "Windows") {
    out <- parallel::mclapply(1:svn_itor, res, mc.cores = 1)
  } else {
    out <- parallel::mclapply(1:svn_itor, res, mc.cores = 1)
  }

  nusvm <- numeric(svn_itor)
  corrv <- numeric(svn_itor)

  for (t in seq_len(svn_itor)) {
    weights <- t(out[[t]]$coefs) %*% out[[t]]$SV
    weights[which(weights < 0)] <- 0
    w <- weights / sum(weights)
    u <- sweep(X, MARGIN = 2, w, "*")
    k <- apply(u, 1, sum)
    nusvm[t] <- sqrt(mean((k - y)^2))
    corrv[t] <- cor(k, y)
  }

  mn <- which.min(nusvm)
  model <- out[[mn]]

  q <- t(model$coefs) %*% model$SV
  q[which(q < 0)] <- 0

  if (!absolute || abs_method == "sig.score") {
    w <- q / sum(q)
  } else if (absolute && abs_method == "no.sumto1") {
    w <- q
  }

  list(w = w, mix_rmse = nusvm[mn], mix_r = corrv[mn])
}

#' Permutation Test for CIBERSORT
#'
#' @description
#' Performs permutation-based sampling to generate an empirical null distribution
#' of correlation coefficients for p-value calculation in CIBERSORT analysis.
#' Randomly samples from the mixture data to create null distributions.
#'
#' @param perm Integer. Number of permutations to perform (\eqn{\ge}100 recommended for
#'   reliable p-value estimation).
#' @param X Matrix or data frame containing signature matrix (predictor variables).
#' @param Y Numeric vector containing the mixture sample expression.
#' @param absolute Logical indicating whether to use absolute space for weights.
#' @param abs_method String specifying the method for absolute space weights:
#'   `"sig.score"` or `"no.sumto1"`.
#' @param seed Integer. Random seed for reproducibility. If NULL (default),
#'   uses current random state.
#'
#' @return List containing:
#' \describe{
#'   \item{dist}{Numeric vector of correlation coefficients from permutations}
#' }
#'
#' @export
#'
#' @examples
#' # Simulate data
#' set.seed(123)
#' X <- matrix(rnorm(100), nrow = 10)
#' rownames(X) <- paste0("Gene", 1:10)
#' colnames(X) <- paste0("Cell", 1:10)
#' Y <- rnorm(10)
#' names(Y) <- paste0("Gene", 1:10)
#'
#' # Run permutation test
#' result <- doPerm(10, X, Y, absolute = FALSE, abs_method = "sig.score")
#' if (!is.null(result)) str(result)
doPerm <- function(perm, X, Y, absolute, abs_method, seed = NULL) {
  if (is.null(perm)) return(NULL)
  if (!is.null(seed)) set.seed(seed)

  if (perm < 1) {
    cli::cli_abort("{.arg perm} must be at least 1, got {perm}")
  }

  Ylist <- as.list(data.matrix(Y))
  dist <- numeric(perm)

  for (itor in seq_len(perm)) {
    yr <- as.numeric(Ylist[sample(length(Ylist), dim(X)[1])])
    yr <- (yr - mean(yr)) / sd(yr)
    result <- CoreAlg(X, yr, absolute, abs_method)
    dist[itor] <- result$mix_r
  }

  list(dist = dist)
}


#' Parallel Permutation Test for CIBERSORT
#'
#' @description
#' Parallel version of doPerm. Performs permutation-based sampling and runs
#' the CoreAlg function iteratively using multiple CPU cores to accelerate
#' computation. This function generates an empirical null distribution of
#' correlation coefficients for p-value calculation in CIBERSORT analysis.
#'
#' @param perm1 Integer. Number of permutations to perform (\eqn{\ge}100 recommended).
#' @param X1 Matrix or data frame. Signature matrix (cell type GEP barcode).
#' @param Y1 Matrix. Mixture file containing gene expression profiles.
#' @param absolute1 Logical. Whether to run in absolute mode (default: FALSE).
#' @param abs_method1 String. Method for absolute mode: 'sig.score' or 'no.sumto1'.
#' @param num_cores1 Integer. Number of CPU cores for parallel computation (default: 2).
#' @param seed Integer. Random seed for reproducibility. If NULL (default),
#'   uses current random state. Set to a specific value (e.g., 123) for
#'   reproducible results across runs.
#'
#' @return List containing:
#' \describe{
#'   \item{dist}{Numeric vector of correlation coefficients from permutations,
#'     representing the empirical null distribution.}
#' }
#'
#' @details
#' This function utilizes the \code{foreach} and \code{doParallel} packages
#' to distribute permutation iterations across multiple cores. It automatically
#' handles cluster setup/teardown via \code{on.exit()} to prevent resource leaks.
#'
#' Note: Windows users may experience slower performance due to socket-based
#' parallelization (PSOCK) versus forking on Unix systems.
#'
#' @seealso \code{\link{doPerm}} for the sequential version, \code{\link{CoreAlg}}, \code{\link{CIBERSORT}}
#'
#' @importFrom parallel makeCluster stopCluster
#' @export
#'
#' @examples
#' # Simulate data
#' set.seed(123)
#' X <- matrix(rnorm(1000), nrow = 100)
#' Y <- matrix(rnorm(500), nrow = 100)
#' rownames(X) <- rownames(Y) <- paste0("Gene", 1:100)
#' colnames(X) <- paste0("Cell", 1:10)
#' colnames(Y) <- paste0("Sample", 1:5)
#'
#' # Run parallel permutation test
#' result <- parallel_doperm(
#'   perm1 = 10, X1 = X, Y1 = Y,
#'   absolute1 = FALSE, abs_method1 = "sig.score", num_cores1 = 1
#' )
#' if (!is.null(result)) str(result$dist)
parallel_doperm <- function(perm1, X1, Y1, absolute1, abs_method1,
                            num_cores1 = 2, seed = NULL) {
  if (is.null(perm1)) return(NULL)
  rlang::check_installed("foreach", reason = "for parallel permutation")
  rlang::check_installed("doParallel", reason = "for parallel backend")

  if (perm1 < 1) {
    cli::cli_abort("{.arg perm1} must be at least 1, got {perm1}")
  }
  if (num_cores1 < 1) {
    cli::cli_abort("{.arg num_cores1} must be at least 1, got {num_cores1}")
  }

  `%dopar%` <- foreach::`%dopar%`

  if (!is.null(seed)) set.seed(seed)
  n <- dim(X1)[1]
  Ylist_len <- length(as.list(data.matrix(Y1)))
  perm_indices <- lapply(seq_len(perm1), function(i) sample(Ylist_len, n))

  cl <- parallel::makeCluster(num_cores1)
  on.exit(
    {
      parallel::stopCluster(cl)
      foreach::registerDoSEQ()
    },
    add = TRUE
  )

  doParallel::registerDoParallel(cl)

  dist <- foreach::foreach(
    itor = seq_len(perm1),
    .combine = rbind,
    .export = c("perm_indices", "Y1", "X1", "absolute1", "abs_method1"),
    .packages = "e1071"
  ) %dopar% {
    Ylist <- as.list(data.matrix(Y1))
    yr <- as.numeric(Ylist[perm_indices[[itor]]])
    yr <- (yr - mean(yr)) / sd(yr)

    result <- CoreAlg(X1, yr, absolute1, abs_method1)
    result$mix_r
  }

  list(dist = dist)
}


#' CIBERSORT Deconvolution Algorithm
#'
#' An analytical tool to estimate cell type abundances in mixed cell populations
#' using gene expression data.

#' @param sig_matrix Cell type GEP barcode matrix: row 1 = sample labels;
#'   column 1 = gene symbols; no missing values; default = LM22.txt download
#'   from CIBERSORT (https://cibersortx.stanford.edu/runcibersort.php)
#' @param mixture_file GEP matrix: row 1 = sample labels; column 1 = gene
#'   symbols; no missing values
#' @param perm Set permutations for statistical analysis
#'   (>=100 permutations recommended).
#' @param QN Quantile normalization of input mixture (default = TRUE)
#' @param absolute Run CIBERSORT in absolute mode (default = FALSE)
#'   - note that cell subsets will be scaled by their absolute levels and
#'     will not be represented as fractions (to derive the default output,
#'     normalize absolute levels such that they sum to 1 for each mixture sample)
#'   - the sum of all cell subsets in each mixture sample will be added to the
#'     output ('Absolute score'). If LM22 is used, this score will capture
#'     total immune content.
#' @param abs_method If absolute is set to TRUE, choose method:
#'   'no.sumto1' or 'sig.score'
#'   - sig.score = for each mixture sample, define S as the median expression
#'     level of all genes in the signature matrix divided by the median
#'     expression level of all genes in the mixture. Multiple cell subset
#'     fractions by S.
#'   - no.sumto1 = remove sum to 1 constraint
#' @param parallel Logical. Enable parallel execution? (default = FALSE)
#' @param num_cores Integer. Number of cores to use when \code{parallel = TRUE}
#'   (default = 2)
#' @param seed Integer. Random seed for reproducible permutation testing.
#'   If \code{NULL} (default), uses current random state. Set to a specific
#'   value (e.g., 123) for reproducible results across runs. Applies to both
#'   parallel and serial permutation.
#' @author Aaron M. Newman, Stanford University (amnewman@stanford.edu)
#' @return A matrix object containing the estimated cibersort-cell fractions,
#'   p-values, correlation coefficients, and RMSE values.
#' @export
#' @import parallel
#' @import dplyr
#' @import tibble
#' @import tidyr
#' @import purrr
#' @import stringr
#' @examples
#' # Simulate data
#' set.seed(123)
#' sim_sig <- matrix(rnorm(100 * 5), 100, 5)
#' rownames(sim_sig) <- paste0("Gene", 1:100)
#' colnames(sim_sig) <- paste0("Cell", 1:5)
#'
#' sim_mixture <- as.data.frame(matrix(
#'   rnorm(100 * 3), 100, 3
#' ))
#' rownames(sim_mixture) <- paste0("Gene", 1:100)
#' colnames(sim_mixture) <- paste0("Sample", 1:3)
#'
#' # Run deconvolution
#' result <- CIBERSORT(
#'   sig_matrix = sim_sig,
#'   mixture_file = sim_mixture,
#'   perm = 10, QN = FALSE, absolute = FALSE,
#'   parallel = FALSE
#' )
#' if (!is.null(result)) head(result)
CIBERSORT <- function(sig_matrix = NULL, mixture_file, perm, QN = TRUE,
                      absolute = FALSE, abs_method = "sig.score",
                      parallel = FALSE, num_cores = 2, seed = NULL) {
  if (is.null(mixture_file)) return(NULL)
  if (is.null(sig_matrix)) {
    sig_matrix <- load_data("lm22")
    if (is.null(sig_matrix)) return(NULL)
  }
  # Input validation
  if (length(intersect(rownames(mixture_file), rownames(sig_matrix))) == 0) {
    stop("None identical gene between eset and reference had been found.
         Check your eset using: intersect(rownames(eset), rownames(reference))")
  }

  if (absolute && abs_method != "no.sumto1" && abs_method != "sig.score") {
    stop("abs_method must be set to either 'sig.score' or 'no.sumto1'")
  }

  # Read in data
  X <- sig_matrix
  Y <- rownames_to_column(mixture_file, var = "symbol")

  # Handle duplicated gene symbols
  dups <- dim(Y)[1] - length(unique(Y[, 1]))
  if (dups > 0) {
    warning(paste(dups, "duplicated gene symbol(s) found in mixture file!"))
    rownames(Y) <- make.names(Y[, 1], unique = TRUE)
  } else {
    rownames(Y) <- Y[, 1]
  }
  Y <- Y[, -1, drop = FALSE]

  X <- data.matrix(X)
  Y <- data.matrix(Y)

  # Order
  X <- X[order(rownames(X)), , drop = FALSE]
  Y <- Y[order(rownames(Y)), , drop = FALSE]

  P <- perm

  # Anti-log if max < 50 in mixture file
  if (max(Y) < 50) {
    Y <- 2^Y
  }

  # Quantile normalization of mixture file
  if (QN == TRUE) {
    tmpc <- colnames(Y)
    tmpr <- rownames(Y)
    rlang::check_installed("preprocessCore")
    Yq <- qn_matrix(Y)
    Y <- matrix(Yq,
      nrow = length(tmpr), ncol = length(tmpc),
      dimnames = list(tmpr, tmpc)
    )
  }

  # Store original mixtures
  Yorig <- Y
  Ymedian <- max(median(Yorig), 1)

  # Intersect genes
  Xgns <- row.names(X)
  Ygns <- row.names(Y)
  YintX <- Ygns %in% Xgns
  Y <- Y[YintX, , drop = FALSE]
  XintY <- Xgns %in% row.names(Y)
  X <- X[XintY, , drop = FALSE]

  # Sanity check
  if (length(dim(Y)) != 2) {
    stop(
      "`mixture_file` became non-2D after preprocessing. ",
      "This usually occurs when it has only one column and subsetting ",
      "dropped dimensions. ",
      "Please update package or ensure `mixture_file` is a matrix with drop = FALSE."
    )
  }

  # Standardize sig matrix
  X <- (X - mean(X)) / sd(as.vector(X))

  # Parallel setup with safety checks
  use_parallel_mixtures <- parallel

  if (parallel) {
    rlang::check_installed("BiocParallel", reason = "for parallel mode")

    # Cap cores to prevent system overload
    max_available <- parallel::detectCores()
    if (is.na(max_available)) max_available <- 2
    safe_cores <- max(1, max_available - 1)

    if (num_cores > safe_cores) {
      warning(
        "Requested ", num_cores, " cores but only ", safe_cores,
        " available/safe. Reducing to ", safe_cores, "."
      )
      num_cores <- safe_cores
    }

    BPPARAM <- BiocParallel::MulticoreParam(workers = num_cores)
  }

  # Empirical null distribution of correlation coefficients
  if (P > 0) {
    if (parallel && exists("parallel_doperm", mode = "function")) {
      # Use parallel permutation (consumes parallel resources)
      kk <- parallel_doperm(
        perm1 = P, X1 = X, Y1 = Y,
        absolute1 = absolute, abs_method1 = abs_method,
        num_cores1 = num_cores, seed = seed
      )
      nulldist <- sort(kk$dist)
      # CRITICAL: Disable parallel for mixtures to avoid nested parallelism
      use_parallel_mixtures <- FALSE
    } else {
      # Serial permutation
      nulldist <- sort(doPerm(P, X, Y, absolute, abs_method, seed = seed)$dist)
    }
  } else {
    nulldist <- NULL
  }

  header <- c("Mixture", colnames(X), "P-value", "Correlation", "RMSE")
  if (absolute) {
    header <- c(header, paste0("Absolute score (", abs_method, ")"))
  }

  # Core job function with explicit parameters (Windows-safe)
  coreJob <- function(itor, Y, X, absolute, abs_method, Ymedian, P, nulldist) {
    y <- Y[, itor, drop = TRUE]
    y <- (y - mean(y)) / sd(y)
    result <- CoreAlg(X, y, absolute, abs_method)
    w <- result$w
    mix_r <- result$mix_r
    mix_rmse <- result$mix_rmse
    if (absolute && abs_method == "sig.score") {
      w <- w * median(Y[, itor]) / Ymedian
    }
    pval <- if (P > 0) {
      1 - which.min(abs(nulldist - mix_r)) / length(nulldist)
    } else {
      9999
    }

    out <- c(colnames(Y)[itor], w, pval, mix_r, mix_rmse)
    if (absolute) out <- c(out, sum(w))
    out
  }

  mixtures <- ncol(Y)

  # Process mixtures
  if (use_parallel_mixtures && mixtures > 1) {
    # Parallel processing (Windows-safe via explicit variable passing)
    output_list <- BiocParallel::bplapply(
      seq_len(mixtures),
      function(i) {
        coreJob(i,
          Y = Y, X = X, absolute = absolute,
          abs_method = abs_method, Ymedian = Ymedian,
          P = P, nulldist = nulldist
        )
      },
      BPPARAM = BPPARAM
    )
    output <- do.call(rbind, output_list)
  } else {
    # Serial processing (safe fallback)
    output_list <- lapply(
      seq_len(mixtures),
      function(i) {
        coreJob(i,
          Y = Y, X = X, absolute = absolute,
          abs_method = abs_method, Ymedian = Ymedian,
          P = P, nulldist = nulldist
        )
      }
    )
    output <- do.call(rbind, output_list)
  }

  # Return matrix object containing all results
  obj <- rbind(header, output)
  obj <- obj[, -1, drop = FALSE]
  obj <- obj[-1, , drop = FALSE]
  obj <- matrix(as.numeric(unlist(obj)), nrow = nrow(obj))
  rownames(obj) <- colnames(Y)
  if (!absolute) {
    colnames(obj) <- c(colnames(X), "P-value", "Correlation", "RMSE")
  } else {
    colnames(obj) <- c(
      colnames(X), "P-value", "Correlation", "RMSE",
      paste0("Absolute score (", abs_method, ")")
    )
  }
  obj
}


############################################################
# FILE: step10_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step10_sensitivity.R


############################################################
# FILE: step10_sensitivity.R
############################################################
# Step 10: Sensitivity analysis - rebuild model WITHOUT MIR644A (reviewer request)
# Rationale: (a) original paper's GRG list used GeneCards "Protein Coding" filter,
#            which excludes miRNAs; our GO/UniProt path lacked that filter.
#            (b) MIR644A cannot be validated cross-platform (no probe on GPL571).
# Decision: rebuild the main model without MIR644A; move MIR644A to discussion.
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(glmnet); library(randomForest); library(pROC) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step10_sensitivity")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
grdegs <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
gset <- setdiff(grdegs$gene, "MIR644A")           # 27 genes
gset <- intersect(gset, rownames(expr))
X <- t(expr[gset, ])
y <- factor(pheno$group, levels = c("Control", "DN"))
cat("27-GRDEG model input:", dim(X), "\n")

# RF (identical settings to original pipeline)
set.seed(2024)
rf <- randomForest(x = X, y = y, ntree = 1000, importance = TRUE)
imp <- importance(rf, type = 2, scale = FALSE)
imp <- data.frame(gene = rownames(imp), MeanDecreaseGini = as.numeric(imp[, 1]))
imp <- imp[order(-imp$MeanDecreaseGini), ]
top10 <- head(imp, 10)$gene
cat("RF top10 (no MIR644A):", paste(top10, collapse = ","), "\n")
write.csv(imp, file.path(outdir, "RF_importance_noMIR644A.csv"), row.names = FALSE)

# LASSO
X10 <- X[, top10]
cvfit <- cv.glmnet(x = X10, y = y, family = "binomial", alpha = 1,
                   nfolds = 10, standardize = TRUE)
cf0 <- as.matrix(coef(cvfit, s = "lambda.min"))
nz <- data.frame(gene = rownames(cf0), coefficient = as.numeric(cf0))
nz <- nz[nz$coefficient != 0 & nz$gene != "(Intercept)", ]
cat("LASSO model genes (no MIR644A, lambda.min =", signif(cvfit$lambda.min, 4), "):\n")
print(nz)
write.csv(nz, file.path(outdir, "LASSO_model_genes_noMIR644A.csv"), row.names = FALSE)
mg <- nz$gene
writeLines(mg, file.path(outdir, "model_genes_noMIR644A.txt"))
rs <- as.numeric(as.matrix(X[, mg]) %*% nz$coefficient[match(mg, nz$gene)])
write.csv(data.frame(gsm = pheno$gsm, group = pheno$group, RiskScore = rs),
          file.path(outdir, "riskscore_noMIR644A_training.csv"), row.names = FALSE)

# training AUC
r1 <- roc(response = y, predictor = rs, quiet = TRUE)
ci1 <- ci.auc(r1)
auc_train <- auc(r1)
cat(sprintf("Training AUC (no MIR644A) = %.4f (%.4f-%.4f)\n", auc_train, ci1[1], ci1[3]))

# external validation with the new model (handle cross-platform missing genes)
auc_rows <- data.frame(dataset = "GSE96804_training", predictor = "RiskScore_noMIR644A",
                       AUC = auc_train, CI_low = ci1[1], CI_high = ci1[3])
for (gse in c("GSE30528", "GSE30529")) {
  ev <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                 check.names = FALSE, stringsAsFactors = FALSE)
  rownames(ev) <- ev$gene; ev$gene <- NULL; ev <- as.matrix(ev)
  phv <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  stopifnot(identical(phv$gsm, colnames(ev)))
  yv <- factor(phv$group, levels = c("Control", "DKD"))
  present <- mg[mg %in% rownames(ev)]
  missing <- setdiff(mg, present)
  cat(gse, "present:", paste(present, collapse = ","), " missing:", paste(missing, collapse = ","), "\n")
  rsv <- as.numeric(t(ev[present, ]) %*% nz$coefficient[match(present, nz$gene)])
  rv <- roc(response = yv, predictor = rsv, quiet = TRUE)
  cv2 <- ci.auc(rv)
  cat(sprintf("%s AUC (no MIR644A) = %.4f (%.4f-%.4f)\n", gse, auc(rv), cv2[1], cv2[3]))
  auc_rows <- rbind(auc_rows, data.frame(dataset = gse, predictor = "RiskScore_noMIR644A",
                                         AUC = auc(rv), CI_low = cv2[1], CI_high = cv2[3]))
}
write.csv(auc_rows, file.path(outdir, "AUC_noMIR644A_model.csv"), row.names = FALSE)

# side-by-side comparison with the original 7-gene model
cmp <- data.frame(
  model = c("Original 7-gene (MIR644A+6)", "Sensitivity 27-GRDEG (no MIR644A)"),
  training_AUC = c(0.9988, auc_train),
  GSE30528_AUC = c(0.7949, auc_rows$AUC[auc_rows$dataset == "GSE30528"]),
  GSE30529_AUC = c(0.7083, auc_rows$AUC[auc_rows$dataset == "GSE30529"]),
  model_genes = c("MIR644A, FN1, SLIT3, LPL, CCL21, SELL, EGF", paste(mg, collapse = ", "))
)
write.csv(cmp, file.path(outdir, "model_comparison_original_vs_noMIR644A.csv"), row.names = FALSE)
cat("\nCOMPARISON:\n"); print(cmp, digits = 3)
cat("\nSTEP10 DONE\n")


############################################################
# FILE: step11_immune_decontam.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 11: Reviewer-response analyses for immune infiltration.
(i) Signature decontamination: remove ALL model genes from the 22 LM22-derived
    signatures, re-run ssGSEA, recompute gene-cell correlations.
(iii) Unified BH (FDR) correction over all 154 gene x cell-type Spearman tests
    (both original and decontaminated), replacing the earlier uncorrected p<0.05.
Extra mechanism check: NeutrophilScore built from classic neutrophil markers
    (LM22 neutrophil signature EXCLUDING SELL) -> correlation with SELL expression.
"""
import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.stats import spearmanr
import gseapy

BASE = "PROJECT_ROOT"
OUT = os.path.join(BASE, "output", "step11_immune_decontam")
os.makedirs(OUT, exist_ok=True)

def bh_adjust(pvals):
    p = np.asarray(pvals, dtype=float)
    n = len(p)
    order = np.argsort(p)
    ranked = p[order] * n / np.arange(1, n + 1)
    q = np.minimum.accumulate(ranked[::-1])[::-1]
    q = np.minimum(q, 1.0)
    out = np.empty(n)
    out[order] = q
    return out

expr = pd.read_csv(os.path.join(BASE, "data", "GSE96804_norm_expr.csv"), index_col=0)
expr.columns = [c.strip() for c in expr.columns]
lm22 = pd.read_csv(os.path.join(BASE, "data", "LM22_signature.csv"), index_col=0)
cell_types = list(lm22.columns)

old_model = [g.strip() for g in open(os.path.join(BASE, "output/step6/model_genes.txt"))]
new_model = [g.strip() for g in open(os.path.join(BASE, "output/step10_sensitivity/model_genes_noMIR644A.txt"))]
all_model_genes = sorted(set(old_model) | set(new_model))
print("model genes union:", all_model_genes)

# how many model genes sit inside the signatures? (contamination audit)
contam = {g: list(lm22.columns[lm22.loc[g] == 1]) for g in all_model_genes if g in lm22.index}
print("contaminating entries (model gene in signature):", contam)

# ---- (i) decontaminated signatures ----
lm22_decon = lm22.drop(index=[g for g in all_model_genes if g in lm22.index])
gene_sets = {ct: sorted(set(lm22_decon.index[lm22_decon[ct] == 1])) for ct in cell_types}
res = gseapy.ssgsea(data=expr, gene_sets=gene_sets, outdir=None,
                    sample_norm_method="rank", min_size=5, no_plot=True)
r2 = res.res2d
nes = r2.pivot(index="Name", columns="Term", values="NES")
nes.index = [str(i).strip() for i in nes.index]
nes = nes.loc[[c for c in expr.columns if c in nes.index]]
print("decontaminated NES:", nes.shape)
nes.to_csv(os.path.join(OUT, "immune_nes_matrix_decontaminated.csv"))

# original NES from step9
nes_orig = pd.read_csv(os.path.join(BASE, "output/step9/immune_nes_matrix.csv"), index_col=0)
nes_orig.index = [str(i).strip() for i in nes_orig.index]

def corr_table(genes, nes_mat):
    rows = []
    for g in genes:
        if g not in expr.index:
            continue
        for ct in nes_mat.columns:
            r, p = spearmanr(expr.loc[g, nes_mat.index].values, nes_mat[ct].values)
            rows.append((g, ct, r, p))
    t = pd.DataFrame(rows, columns=["gene", "cell_type", "spearman_r", "p_value"])
    t["q_BH"] = bh_adjust(t["p_value"].values)
    return t

def clean_nes(nes_mat):
    m = nes_mat.copy()
    m.columns = [str(c).strip() for c in m.columns]
    return m

# A) old model genes x decontaminated NES  (direct fix of Supplementary Table S3)
old_decon = clean_nes(corr_table([g for g in old_model if g in expr.index], nes))
old_decon.to_csv(os.path.join(OUT, "corr_pval_BH_oldmodel_decontaminated.csv"), index=False)

# B) new (final) model genes x decontaminated NES
new_decon = clean_nes(corr_table([g for g in new_model if g in expr.index], nes))
new_decon.to_csv(os.path.join(OUT, "corr_pval_BH_newmodel_decontaminated.csv"), index=False)

# original tables recomputed with BH for consistency
old_orig = clean_nes(corr_table([g for g in old_model if g in expr.index], nes_orig))
old_orig.to_csv(os.path.join(OUT, "corr_pval_BH_oldmodel_originalNES.csv"), index=False)

# ---- key comparisons ----
def get(t, g, ct):
    row = t[(t.gene == g) & (t.cell_type == ct)]
    return row.iloc[0] if len(row) else None

cmp_rows = []
for g in sorted(set(old_model) | set(new_model)):
    if g not in expr.index:
        continue
    for ct in cell_types:
        o = get(old_orig, g, ct)
        d = get(old_decon, g, ct)
        d2 = get(new_decon, g, ct)
        if o is not None and d is not None:
            cmp_rows.append(dict(gene=g, cell_type=ct,
                                 r_original=o.spearman_r, p_original=o.p_value, q_original=o.q_BH,
                                 r_decontaminated=d.spearman_r, p_decontaminated=d.p_value, q_decontaminated=d.q_BH,
                                 r_newmodel_decontaminated=(d2.spearman_r if d2 is not None else np.nan),
                                 q_newmodel_decontaminated=(d2.q_BH if d2 is not None else np.nan)))
cmp = pd.DataFrame(cmp_rows)
cmp.to_csv(os.path.join(OUT, "original_vs_decontaminated_key_pairs.csv"), index=False)

sn_orig = get(old_orig, "SELL", "Neutrophils")
sn_decon = get(old_decon, "SELL", "Neutrophils")
sb_orig = get(old_orig, "SELL", "B cells naive")
sb_decon = get(old_decon, "SELL", "B cells naive")
print("\nSELL-Neutrophils: original r=%.3f (q=%.2e) -> decontaminated r=%.3f (q=%.2e)"
      % (sn_orig.spearman_r, sn_orig.q_BH, sn_decon.spearman_r, sn_decon.q_BH))
print("SELL-Bnaive(signature-member): original r=%.3f -> decontaminated r=%.3f"
      % (sb_orig.spearman_r, sb_decon.spearman_r))

# ---- mechanism check: neutrophil biology vs SELL (no signature inference) ----
neut_genes = sorted(set(lm22.index[lm22["Neutrophils"] == 1]))
neut_genes_noSELL = [g for g in neut_genes if g != "SELL" and g in expr.index]
zs = expr.loc[neut_genes_noSELL].apply(lambda r: (r - r.mean()) / (r.std() + 1e-12), axis=1)
neut_score = zs.mean(axis=0)
r_bio, p_bio = spearmanr(expr.loc["SELL", neut_score.index].values, neut_score.values)
rsell, psell = spearmanr(expr.loc["SELL", nes.index].values, nes_orig.loc[nes.index, "Neutrophils"].values)
print("NeutrophilScore (markers excl. SELL, n=%d) vs SELL expression: r=%.3f p=%.2e"
      % (len(neut_genes_noSELL), r_bio, p_bio))
pd.DataFrame(dict(neutrophil_markers_excluding_SELL=[neut_genes_noSELL[0]] + [""] * (len(neut_genes_noSELL) - 1),
                  all=list(neut_genes_noSELL))).to_csv(
    os.path.join(OUT, "neutrophil_markers_used.csv"), index=False)
with open(os.path.join(OUT, "mechanism_check_summary.txt"), "w") as f:
    f.write("SELL-NeutrophilsNES original: r=%.4f p=%.3e\n" % (rsell, psell))
    f.write("SELL-NeutrophilScore(marker-mean z, excl SELL): r=%.4f p=%.3e\n" % (r_bio, p_bio))
    f.write("SELL-BnaiveNES original r=%.4f -> decontaminated r=%.4f (signature-member auto-correlation)\n"
            % (sb_orig.spearman_r, sb_decon.spearman_r))
    f.write("SELL-NeutrophilsNES original r=%.4f -> decontaminated r=%.4f\n"
            % (sn_orig.spearman_r, sn_decon.spearman_r))

# ---- significance counts under BH (unified statistics) ----
sig_counts = pd.DataFrame([
    dict(table="old model x original NES (Suppl Table S3 fix)",
         n_tests=len(old_orig), n_sig_p005=int((old_orig.p_value < 0.05).sum()),
         n_sig_q005=int((old_orig.q_BH < 0.05).sum())),
    dict(table="old model x decontaminated NES",
         n_tests=len(old_decon), n_sig_p005=int((old_decon.p_value < 0.05).sum()),
         n_sig_q005=int((old_decon.q_BH < 0.05).sum())),
    dict(table="new model x decontaminated NES (final)",
         n_tests=len(new_decon), n_sig_p005=int((new_decon.p_value < 0.05).sum()),
         n_sig_q005=int((new_decon.q_BH < 0.05).sum())),
])
sig_counts.to_csv(os.path.join(OUT, "significance_counts_BH_unified.csv"), index=False)
print("\n", sig_counts)

# ---- comparison figure ----
fig, axes = plt.subplots(1, 2, figsize=(12, 4.6))
sub = cmp[(cmp.gene == "SELL")]
x = np.arange(len(sub))
axes[0].bar(x - 0.2, sub.r_original, width=0.4, label="original NES", color="#B2182B", alpha=0.85)
axes[0].bar(x + 0.2, sub.r_decontaminated, width=0.4, label="decontaminated NES", color="#2166AC", alpha=0.85)
axes[0].set_xticks(x)
axes[0].set_xticklabels([c.replace(" cells", "").replace("T cells ", "T-") for c in sub.cell_type],
                        rotation=90, fontsize=6)
axes[0].set_ylabel("Spearman r with SELL")
axes[0].set_title("SELL vs 22 immune cells: original vs decontaminated")
axes[0].legend(fontsize=8)
axes[0].axhline(0, color="grey", lw=0.6)
top = new_decon.reindex(new_decon.r_spearman_r.abs().sort_values().index
                        if hasattr(new_decon, "r_spearman_r") else new_decon.spearman_r.abs().sort_values().index)
top_sig = new_decon[new_decon.q_BH < 0.05].sort_values("spearman_r")
y = np.arange(len(top_sig))
axes[1].barh(y, top_sig.spearman_r, color=["#2166AC" if v < 0 else "#B2182B" for v in top_sig.spearman_r])
axes[1].set_yticks(y)
axes[1].set_yticklabels([f"{g} | {c}" for g, c in zip(top_sig.gene, top_sig.cell_type)], fontsize=6)
axes[1].set_xlabel("Spearman r (q_BH < 0.05 only)")
axes[1].set_title("Final model genes vs immune cells (BH-significant)")
plt.tight_layout()
plt.savefig(os.path.join(OUT, "fig_decontam_comparison.svg"), bbox_inches="tight")
plt.close()
print("\nSTEP11 DONE")


############################################################
# FILE: step13_formal_rerun.R
############################################################
# Step 13: Formal rerun with 27 protein-coding GRDEGs (lambda.min + lambda.1se)
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(glmnet); library(randomForest); library(pROC) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step13_formal_rerun")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
grd <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
gset <- setdiff(grd$gene, "MIR644A")            # 27 protein-coding GRDEGs
X <- t(expr[gset, ])
y <- factor(pheno$group, levels = c("Control", "DN"))

# ---- univariate logistic (Suppl Table S2 v2) ----
uni <- do.call(rbind, lapply(gset, function(g) {
  d <- data.frame(y = y, x = X[, g])
  fit <- suppressWarnings(glm(y ~ x, family = binomial, data = d))
  s <- summary(fit)$coefficients
  data.frame(gene = g, OR = exp(s["x", 1]), CI_low = exp(s["x", 1] - 1.96 * s["x", 2]),
             CI_high = exp(s["x", 1] + 1.96 * s["x", 2]), p = s["x", 4],
             logFC = grd$logFC[match(g, grd$gene)], stringsAsFactors = FALSE)
}))
uni <- uni[order(uni$p), ]
write.csv(uni, file.path(outdir, "S2_univariate_logistic_27GRDEG.csv"), row.names = FALSE)
cat("univariate significant (p<0.05):", sum(uni$p < 0.05), "/ 27\n")

# ---- RF top10 ----
set.seed(2024)
rf <- randomForest(x = X, y = y, ntree = 1000, importance = TRUE)
imp <- importance(rf, type = 2, scale = FALSE)
imp <- data.frame(gene = rownames(imp), MeanDecreaseGini = as.numeric(imp[, 1]))
imp <- imp[order(-imp$MeanDecreaseGini), ]
top10 <- head(imp, 10)$gene
write.csv(imp, file.path(outdir, "RF_importance_27GRDEG.csv"), row.names = FALSE)
cat("RF top10:", paste(top10, collapse = ","), "\n")

# ---- LASSO at lambda.min AND lambda.1se ----
X10 <- X[, top10]
cvfit <- cv.glmnet(x = X10, y = y, family = "binomial", alpha = 1,
                   nfolds = 10, standardize = TRUE)
models <- list()
for (lmd in c("lambda.min", "lambda.1se")) {
  cf0 <- as.matrix(coef(cvfit, s = lmd))
  nz <- data.frame(gene = rownames(cf0), coefficient = as.numeric(cf0))
  nz <- nz[nz$coefficient != 0 & nz$gene != "(Intercept)", ]
  w <- rep(0, ncol(X10)); names(w) <- colnames(X10)
  w[nz$gene] <- nz$coefficient
  rs <- as.numeric(X10 %*% w)
  r <- roc(response = y, predictor = rs, quiet = TRUE)
  ci <- ci.auc(r)
  models[[lmd]] <- list(genes = nz, rs = rs, auc = auc(r), ci = ci)
  cat(sprintf("%s: lambda=%.5f  genes: %s\n   training AUC=%.4f (%.4f-%.4f)\n",
              lmd, cvfit[[lmd]], paste(nz$gene, collapse = ","), auc(r), ci[1], ci[3]))
  write.csv(nz, file.path(outdir, paste0("LASSO_genes_", lmd, ".csv")), row.names = FALSE)
  write.csv(data.frame(gsm = pheno$gsm, group = pheno$group, RiskScore = rs),
            file.path(outdir, paste0("riskscore_training_", lmd, ".csv")), row.names = FALSE)
}
write.csv(data.frame(lambda_type = c("lambda.min", "lambda.1se"),
                     lambda = c(cvfit$lambda.min, cvfit$lambda.1se),
                     AUC = c(models[["lambda.min"]]$auc, models[["lambda.1se"]]$auc)),
          file.path(outdir, "LASSO_min_vs_1se_trainingAUC.csv"), row.names = FALSE)

# main model = lambda.min (7 genes); save formula
mg <- models[["lambda.min"]]$genes
w <- mg$coefficient; names(w) <- mg$gene
formula_txt <- paste0("RiskScore = ", paste0(sprintf("(%.6f x Exp(%s))", w, mg$gene), collapse = " + "))
writeLines(c(formula_txt, paste(mg$gene, collapse = ",")),
           file.path(outdir, "model_formula_and_genes.txt"))
cat("\n", formula_txt, "\n")
saveRDS(cvfit, file.path(outdir, "cv_glmnet.rds"))
cat("STEP13 DONE\n")


############################################################
# FILE: step13_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step13_formal_rerun.R


############################################################
# FILE: step14_robustness.R
############################################################
# Step 14: Model robustness - bootstrap stability + MWU check + repeated 5-fold CV
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(glmnet); library(randomForest); library(pROC) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step14_robustness")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
grd <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
gset <- setdiff(grd$gene, "MIR644A")
Xall <- t(expr[gset, ])
y <- factor(pheno$group, levels = c("Control", "DN"))
n <- length(y)
top10 <- readLines(file.path(base, "output", "step13_formal_rerun", "model_formula_and_genes.txt"))[2]
top10 <- strsplit(top10, ",")[[1]]           # lambda.min 7 genes = formal model
# for bootstrap we start from the 27 genes (full re-selection each draw)

## ---- (A) bootstrap 1000: selection frequency (RF top10 -> LASSO lambda.min) ----
set.seed(2024)
B <- 1000
freq_min <- setNames(numeric(length(gset)), gset)
freq_1se <- setNames(numeric(length(gset)), gset)
coef_sum <- setNames(numeric(length(gset)), gset)
for (b in seq_len(B)) {
  idx <- sample(seq_len(n), replace = TRUE)
  Xb <- Xall[idx, , drop = FALSE]; yb <- y[idx]
  if (length(unique(yb)) < 2) next
  rf <- try(randomForest(x = Xb, y = yb, ntree = 500, importance = FALSE), silent = TRUE)
  if (inherits(rf, "try-error")) next
  im <- importance(rf, type = 2, scale = FALSE)
  im <- im[order(-im[, 1]), , drop = FALSE]
  cand <- head(rownames(im), 10)
  cvb <- try(cv.glmnet(x = Xb[, cand, drop = FALSE], y = yb, family = "binomial",
                       alpha = 1, nfolds = 5, standardize = TRUE), silent = TRUE)
  if (inherits(cvb, "try-error")) next
  for (lmd in c("lambda.min", "lambda.1se")) {
    cf0 <- as.numeric(coef(cvb, s = lmd))
    nm <- rownames(coef(cvb, s = lmd))
    sel <- nm[cf0 != 0 & nm != "(Intercept)"]
    if (lmd == "lambda.min") freq_min[sel] <- freq_min[sel] + 1 else freq_1se[sel] <- freq_1se[sel] + 1
    if (lmd == "lambda.min") coef_sum[sel] <- coef_sum[sel] + cf0[match(sel, nm)]
  }
  if (b %% 200 == 0) cat("bootstrap", b, "done\n")
}
stab <- data.frame(gene = gset,
                   selected_lambda_min = freq_min, freq_lambda_min = freq_min / B,
                   selected_lambda_1se = freq_1se, freq_lambda_1se = freq_1se / B)
stab <- stab[order(-stab$freq_lambda_min), ]
write.csv(stab, file.path(outdir, "bootstrap_selection_frequency_1000.csv"), row.names = FALSE)
cat("\nBootstrap selection frequency (lambda.min, top):\n")
print(head(stab[, c("gene", "freq_lambda_min", "freq_lambda_1se")], 12))

## ---- (B) Mann-Whitney U check on validation RiskScores (old vs new model) ----
mw <- data.frame()
load_val <- function(gse) {
  ev <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                 check.names = FALSE, stringsAsFactors = FALSE)
  rownames(ev) <- ev$gene; ev$gene <- NULL; ev <- as.matrix(ev)
  phv <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  list(ev = ev, phv = phv)
}
old7 <- c("MIR644A", "FN1", "SLIT3", "LPL", "CCL21", "SELL", "EGF")
cf_old <- c(0.93178726, 0.05918551, 1.05350404, -0.96221260, 0.21216499, -1.09470869, -0.48869261)
names(cf_old) <- old7
newtxt <- readLines(file.path(base, "output", "step13_formal_rerun", "model_formula_and_genes.txt"))[2]
new7 <- strsplit(newtxt, ",")[[1]]
cfnew <- read.csv(file.path(base, "output", "step13_formal_rerun", "LASSO_genes_lambda.min.csv"), stringsAsFactors = FALSE)
w_new <- cfnew$coefficient; names(w_new) <- cfnew$gene
for (gse in c("GSE30528", "GSE30529")) {
  d <- load_val(gse); ev <- d$ev; phv <- d$phv
  yv <- factor(phv$group, levels = c("Control", "DKD"))
  present <- old7[old7 %in% rownames(ev)]
  w2 <- cf_old[present]
  rs_old <- as.numeric(t(ev[present, ]) %*% w2)
  rs_new <- as.numeric(t(ev[new7, ]) %*% w_new[new7])
  mt_old <- wilcox.test(rs_old ~ yv, exact = FALSE)
  mt_new <- wilcox.test(rs_new ~ yv, exact = FALSE)
  auc_old <- auc(roc(yv, rs_old, quiet = TRUE)); auc_new <- auc(roc(yv, rs_new, quiet = TRUE))
  mw <- rbind(mw, data.frame(dataset = gse, model = "old7_withMIR644A",
                             AUC = auc_old, MWU_p = mt_old$p.value))
  mw <- rbind(mw, data.frame(dataset = gse, model = "new7_noMIR644A",
                             AUC = auc_new, MWU_p = mt_new$p.value))
  cat(sprintf("%s old AUC=%.4f MWU p=%.4f | new AUC=%.4f MWU p=%.4f\n",
              gse, auc_old, mt_old$p.value, auc_new, mt_new$p.value))
}
write.csv(mw, file.path(outdir, "MWU_validation_check.csv"), row.names = FALSE)

## ---- (C) repeated 5-fold CV (100 repeats), feature re-selection inside each fold ----
set.seed(42)
R <- 100; K <- 5
oof_auc <- numeric(R)
for (r in seq_len(R)) {
  folds <- sample(rep(seq_len(K), length.out = n))
  pred <- numeric(n)
  for (k in seq_len(K)) {
    tr <- which(folds != k); te <- which(folds == k)
    if (length(unique(y[tr])) < 2) next
    rf <- try(randomForest(x = Xall[tr, ], y = y[tr], ntree = 300, importance = FALSE), silent = TRUE)
    if (inherits(rf, "try-error")) next
    im <- importance(rf, type = 2, scale = FALSE)
    im <- im[order(-im[, 1]), , drop = FALSE]
    cand <- head(rownames(im), 10)
    cvk <- try(cv.glmnet(x = Xall[tr, cand, drop = FALSE], y = y[tr], family = "binomial",
                         alpha = 1, nfolds = 4, standardize = TRUE), silent = TRUE)
    if (inherits(cvk, "try-error")) next
    pte <- as.numeric(predict(cvk, newx = Xall[te, cand, drop = FALSE], s = "lambda.min",
                              type = "response"))
    pred[te] <- pte
  }
  if (length(unique(pred)) > 1) oof_auc[r] <- as.numeric(auc(roc(y, pred, quiet = TRUE)))
  if (r %% 20 == 0) cat("CV repeat", r, "OOF AUC =", round(oof_auc[r], 3), "\n")
}
cat(sprintf("\nRepeated 5-fold CV OOF AUC: median=%.4f  IQR=%.4f-%.4f  mean=%.4f  sd=%.4f\n",
            median(oof_auc), quantile(oof_auc, .25), quantile(oof_auc, .75),
            mean(oof_auc), sd(oof_auc)))
write.csv(data.frame(repeat_id = seq_len(R), oof_auc = oof_auc),
          file.path(outdir, "repeated_cv_oof_auc.csv"), row.names = FALSE)
writeLines(sprintf("median=%.4f; IQR=%.4f-%.4f; mean=%.4f; sd=%.4f; n_repeats=%d",
                   median(oof_auc), quantile(oof_auc, .25), quantile(oof_auc, .75),
                   mean(oof_auc), sd(oof_auc), sum(oof_auc > 0)),
           file.path(outdir, "repeated_cv_summary.txt"))

svg(file.path(outdir, "fig_oof_auc_distribution.svg"), width = 7, height = 5)
hist(oof_auc, breaks = 20, col = "#4393C3", border = "white",
     main = "Repeated 5-fold CV: out-of-fold AUC (100 repeats)",
     xlab = "Out-of-fold AUC", cex.main = 0.95)
abline(v = median(oof_auc), lty = 2, col = "#B2182B")
text(median(oof_auc), par("usr")[4] * 0.95, sprintf("median=%.3f", median(oof_auc)),
     pos = 4, cex = 0.9, col = "#B2182B")
dev.off()
cat("STEP14 DONE\n")


############################################################
# FILE: step14_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step14_robustness.R


############################################################
# FILE: step15_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step15_scale_sensitivity.R


############################################################
# FILE: step15_scale_sensitivity.R
############################################################
# Step 15: Cross-platform coefficient-scale sensitivity analysis
# beta were fitted on GPL17586 log2 RMA scale; validation sets are median-centered.
# Sensitivity: recompute RiskScore after WITHIN-COHORT z-scoring of the 7 genes
# (removes scale/offset mismatch; sign of coefficients preserved), compare AUC.
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(pROC) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step15_scale_sensitivity")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

cfnew <- read.csv(file.path(base, "output", "step13_formal_rerun", "LASSO_genes_lambda.min.csv"), stringsAsFactors = FALSE)
w <- cfnew$coefficient; names(w) <- cfnew$gene
rows <- data.frame()
for (gse in c("GSE30528", "GSE30529")) {
  ev <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                 check.names = FALSE, stringsAsFactors = FALSE)
  rownames(ev) <- ev$gene; ev$gene <- NULL; ev <- as.matrix(ev)
  phv <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  yv <- factor(phv$group, levels = c("Control", "DKD"))
  Xg <- t(ev[names(w), ])
  rs_raw <- as.numeric(Xg %*% w)
  Xz <- scale(Xg)                                    # within-cohort z-score
  rs_z <- as.numeric(Xz %*% w)
  Xr <- apply(Xg, 2, rank)                           # within-cohort rank
  rs_r <- as.numeric(Xr %*% w)
  for (m in c("raw(centered-series)", "within-cohort z-score", "within-cohort rank")) {
    rs <- switch(m, "raw(centered-series)" = rs_raw,
                 "within-cohort z-score" = rs_z, "within-cohort rank" = rs_r)
    r <- roc(yv, rs, quiet = TRUE)
    rows <- rbind(rows, data.frame(dataset = gse, scale = m, AUC = as.numeric(auc(r))))
  }
  cat(gse, "raw:", round(auc(roc(yv, rs_raw, quiet = TRUE)), 4),
      " z-score:", round(auc(roc(yv, rs_z, quiet = TRUE)), 4),
      " rank:", round(auc(roc(yv, rs_r, quiet = TRUE)), 4), "\n")
}
write.csv(rows, file.path(outdir, "scale_sensitivity_AUC.csv"), row.names = FALSE)
cat("\nSTEP15 DONE\n")


############################################################
# FILE: step16_downstream.R
############################################################
# Step 16: Downstream recalcs for the 27-GRDEG / 7-gene final model
# (GO/KEGG enrichment, GOSemSim friends, single-gene ROCs incl. FBN1, gene correlation matrix)
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(clusterProfiler); library(org.Hs.eg.db); library(GOSemSim); library(pROC)
  library(ggplot2); library(pheatmap); library(reshape2)
})
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step16_downstream_v2")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

grd <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
genes27 <- setdiff(grd$gene, "MIR644A")
mg <- strsplit(readLines(file.path(base, "output", "step13_formal_rerun", "model_formula_and_genes.txt"))[2], ",")[[1]]

## ---- GO/KEGG enrichment on 27 genes ----
eg <- bitr(genes27, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
for (ont in c("BP", "CC", "MF")) {
  e <- enrichGO(gene = eg$ENTREZID, OrgDb = org.Hs.eg.db, keyType = "ENTREZID",
                ont = ont, pvalueCutoff = 0.05, qvalueCutoff = 0.25,
                pAdjustMethod = "BH", readable = TRUE)
  n <- if (is.null(e)) 0 else nrow(as.data.frame(e))
  cat("GO", ont, "terms:", n, "\n")
  if (n > 0) write.csv(as.data.frame(e), file.path(outdir, paste0("GO_", ont, "_27GRDEG.csv")), row.names = FALSE)
}
link <- read.delim(file.path(base, "data", "kegg_hsa_link.tsv"), header = FALSE, col.names = c("path", "gene"))
link$pathway <- sub("path:", "", link$path); link$entrez <- sub("hsa:", "", link$gene)
pw <- read.delim(file.path(base, "data", "kegg_hsa_pathways.tsv"), header = FALSE)
names(pw) <- c("pathway", "name"); pw$name <- sub(" - Homo sapiens.*$", "", pw$name)
kegg <- enricher(gene = eg$ENTREZID, TERM2GENE = unique(link[, c("pathway", "entrez")]),
                 TERM2NAME = pw[, c("pathway", "name")], pvalueCutoff = 0.05, qvalueCutoff = 0.25)
nK <- if (is.null(kegg) || nrow(as.data.frame(kegg)) == 0) 0 else nrow(as.data.frame(kegg))
cat("KEGG terms:", nK, "\n")
if (nK > 0) {
  kegg <- setReadable(kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
  write.csv(as.data.frame(kegg), file.path(outdir, "KEGG_27GRDEG.csv"), row.names = FALSE)
}

## ---- GOSemSim friends for the final 7-gene model ----
eg7 <- bitr(mg, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
semBP <- godata(org.Hs.eg.db, ont = "BP")
sim <- try(mgeneSim(genes = eg7$ENTREZID, semData = semBP, measure = "Wang"), silent = TRUE)
if (!inherits(sim, "try-error") && !is.null(sim)) {
  sym <- eg7$SYMBOL[match(colnames(sim), eg7$ENTREZID)]
  colnames(sim) <- sym; rownames(sim) <- sym
  keep <- rowSums(!is.na(sim)) > 0 & colSums(!is.na(sim)) > 0
  sim <- sim[keep, keep]
  svg(file.path(outdir, "fig_friends_final7.svg"), width = 7, height = 6)
  pheatmap(sim, display_numbers = TRUE, number_format = "%.2f", number_size = 7,
           color = colorRampPalette(c("#FFF5F0", "#67000D"))(100),
           main = "GO(BP) semantic similarity - final 7-gene model")
  dev.off()
  avg <- sort(rowMeans(sim, na.rm = TRUE), decreasing = TRUE)
  write.csv(data.frame(gene = names(avg), mean_similarity = as.numeric(avg)),
            file.path(outdir, "friends_mean_similarity_final7.csv"), row.names = FALSE)
  cat("friends (final 7):", paste(sprintf("%s=%.3f", names(avg), avg), collapse = " "), "\n")
}

## ---- single-gene ROCs (training + both validation sets) ----
load_expr <- function(gse) {
  e <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                check.names = FALSE, stringsAsFactors = FALSE)
  rownames(e) <- e$gene; e$gene <- NULL; as.matrix(e)
}
rows <- data.frame()
tr <- load_expr("GSE96804"); ytr <- factor(read.csv(file.path(base, "data", "GSE96804_pheno.csv"))$group,
                                           levels = c("Control", "DN"))
for (g in mg) {
  r <- roc(ytr, tr[g, ], quiet = TRUE); c2 <- ci.auc(r)
  rows <- rbind(rows, data.frame(dataset = "GSE96804_training", gene = g,
                                 AUC = as.numeric(auc(r)), CI_low = c2[1], CI_high = c2[3]))
}
for (gse in c("GSE30528", "GSE30529")) {
  ev <- load_expr(gse)
  yv <- factor(read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")))$group,
               levels = c("Control", "DKD"))
  for (g in mg[mg %in% rownames(ev)]) {
    r <- roc(yv, ev[g, ], quiet = TRUE); c2 <- ci.auc(r)
    rows <- rbind(rows, data.frame(dataset = gse, gene = g,
                                   AUC = as.numeric(auc(r)), CI_low = c2[1], CI_high = c2[3]))
  }
}
write.csv(rows, file.path(outdir, "single_gene_AUC_final7.csv"), row.names = FALSE)
print(rows, digits = 3)

## ---- gene-gene correlation matrix (final 7, training set) ----
Xtr <- t(tr[mg, ])
cm <- cor(Xtr)
cmel <- melt(cm); colnames(cmel) <- c("g1", "g2", "r")
write.csv(cmel, file.path(outdir, "gene_correlation_final7.csv"), row.names = FALSE)
pairs <- cmel[cmel$g1 != cmel$g2 & !duplicated(t(apply(cmel[cmel$g1 != cmel$g2, c("g1","g2")], 1, sort))), ]
pairs <- pairs[order(-abs(pairs$r)), ]
cat("strongest gene-gene pairs:\n"); print(head(pairs, 5))
cat("\nSTEP16 DONE\n")


############################################################
# FILE: step16_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step16_downstream.R


############################################################
# FILE: step16b_pairs.R
############################################################
# Step 16b: fix gene-gene pair dedup and print strongest pairs
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(reshape2) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step16_downstream_v2")
mg <- strsplit(readLines(file.path(base, "output", "step13_formal_rerun", "model_formula_and_genes.txt"))[2], ",")[[1]]
e <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
              check.names = FALSE, stringsAsFactors = FALSE)
rownames(e) <- e$gene; e$gene <- NULL
cm <- cor(t(e[mg, ]))
cmel <- melt(cm); colnames(cmel) <- c("g1", "g2", "r")
write.csv(cmel, file.path(outdir, "gene_correlation_final7.csv"), row.names = FALSE)
sub <- cmel[cmel$g1 != cmel$g2, ]
key <- apply(sub[, c("g1", "g2")], 1, function(v) paste(sort(v), collapse = "|"))
pairs <- sub[!duplicated(key), ]
pairs <- pairs[order(-abs(pairs$r)), ]
cat("strongest gene-gene pairs (final 7, training set):\n")
print(head(pairs, 8))
write.csv(pairs, file.path(outdir, "gene_pairs_sorted_final7.csv"), row.names = FALSE)
cat("STEP16B DONE\n")


############################################################
# FILE: step16b_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step16b_pairs.R


############################################################
# FILE: step17_cibersort_svr.R
############################################################
# Step 17: TRUE CIBERSORT nu-SVR deconvolution (official LM22 expression matrix,
# obtained from mdozmorov/Immuno_notes mirror) - answers reviewer (ii) properly.
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(e1071); library(preprocessCore); library(tibble) })
source("PROJECT_ROOT/scripts/CIBERSORT_patched.R")
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step17_cibersort_svr")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL
mix <- 2^as.matrix(expr) - 1                 # back to linear scale
mix <- as.data.frame(mix)
lm22 <- read.delim(file.path(base, "data", "LM22_official.txt"),
                   row.names = 1, check.names = FALSE)
lm22 <- as.matrix(lm22)
cat("mixture:", dim(mix), " LM22:", dim(lm22), "\n")
cat("SELL in LM22 (official expression): row present:", "SELL" %in% rownames(lm22), "\n")
sel <- lm22["SELL", ]
cat("official LM22: top cell types by SELL expression:", paste(names(sort(sel, decreasing = TRUE))[1:3], collapse = ", "), "\n")

# ---- (a) SVR with official LM22 ----
set.seed(2024)
prop <- CIBERSORT(sig_matrix = lm22, mixture_file = mix, perm = 0, QN = TRUE)
stopifnot(!is.null(prop))
prop <- as.data.frame(prop)
prop <- prop[, grep("P.value|Corr|RMSE", colnames(prop), invert = TRUE)]
cat("deconvolved proportions:", dim(prop), "\n")
write.csv(prop, file.path(outdir, "cibersort_proportions_official_LM22.csv"))

# ---- (b) decontaminated LM22 (remove SELL row) -> re-run ----
lm22_d <- lm22[!rownames(lm22) %in% c("SELL"), ]
prop_d <- CIBERSORT(sig_matrix = lm22_d, mixture_file = mix, perm = 0, QN = TRUE)
prop_d <- as.data.frame(prop_d)[, colnames(prop)]
write.csv(prop_d, file.path(outdir, "cibersort_proportions_LM22_noSELL.csv"))

# ---- (c) correlations with SELL expression ----
sell <- expr["SELL", colnames(prop)]
neut_orig <- prop[, "Neutrophils"]
neut_d <- prop_d[, "Neutrophils"]
r_o <- cor.test(sell, neut_orig, method = "spearman")
r_d <- cor.test(sell, neut_d, method = "spearman")
cat(sprintf("SVR(official LM22): SELL-Neutrophil r=%.3f p=%.2e\n", r_o$estimate, r_o$p.value))
cat(sprintf("SVR(LM22 no SELL):  SELL-Neutrophil r=%.3f p=%.2e\n", r_d$estimate, r_d$p.value))
mg <- strsplit(readLines(file.path(base, "output", "step13_formal_rerun", "model_formula_and_genes.txt"))[2], ",")[[1]]
res <- data.frame()
for (g in mg) {
  gexp <- expr[g, colnames(prop)]
  ro <- cor.test(gexp, neut_orig, method = "spearman")
  rd <- cor.test(gexp, neut_d, method = "spearman")
  in_sig <- g %in% rownames(lm22)
  if (in_sig) top_ct <- names(sort(lm22[g, ], decreasing = TRUE))[1] else top_ct <- "not-in-LM22"
  res <- rbind(res, data.frame(gene = g, in_official_LM22 = in_sig, top_LM22_cell = top_ct,
                               r_Neutrophil_origLM22 = ro$estimate, p_orig = ro$p.value,
                               r_Neutrophil_noSELL_LM22 = rd$estimate, p_noSELL = rd$p.value))
}
print(res, digits = 3)
write.csv(res, file.path(outdir, "gene_neutrophil_SVR_corr.csv"), row.names = FALSE)
writeLines(sprintf("SELL-Neutrophil SVR: official LM22 r=%.4f p=%.3e | LM22-without-SELL r=%.4f p=%.3e\ncontamination audit: SELL in official LM22? %s ; top cell by SELL expr: %s",
                   r_o$estimate, r_o$p.value, r_d$estimate, r_d$p.value,
                   "SELL" %in% rownames(lm22),
                   names(sort(lm22["SELL", ], decreasing = TRUE))[1]),
           file.path(outdir, "svr_sell_neutrophil_summary.txt"))

svg(file.path(outdir, "fig_svr_sell_vs_neutrophil.svg"), width = 9, height = 4.5)
par(mfrow = c(1, 2), cex.main = 0.9)
plot(neut_orig, sell, xlab = "Neutrophil fraction (official LM22)",
     ylab = "SELL expression (log2)", pch = 19, col = ifelse(read.csv(file.path(base, "data", "GSE96804_pheno.csv"))$group == "DN", "#E69F00", "#56B4E9"),
     main = sprintf("SVR official LM22: r=%.3f", r_o$estimate))
legend("topleft", c("DN", "Control"), fill = c("#E69F00", "#56B4E9"), cex = 0.8)
plot(neut_d, sell, xlab = "Neutrophil fraction (LM22 without SELL)",
     ylab = "SELL expression (log2)", pch = 19,
     col = ifelse(read.csv(file.path(base, "data", "GSE96804_pheno.csv"))$group == "DN", "#E69F00", "#56B4E9"),
     main = sprintf("SVR no-SELL LM22: r=%.3f", r_d$estimate))
dev.off()
cat("STEP17 DONE\n")


############################################################
# FILE: step17_immune_reinforce.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 17: Immune reinforcement.
(a) 231 pairwise cell-cell correlations with BH FDR;
(b) global-factor audit: PCA on NES, PC1 share, partial correlations controlling PC1;
(c) GO-derived NeutrophilScore (source: GO neutrophil terms, independent of LM22
    construction) with overlap audit vs LM22 neutrophil signature, SELL correlation
    reproduced in all three cohorts;
(d) one more MCP-counter/LM22 mirror attempt (conditional conclusion).
"""
import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.stats import spearmanr
import gseapy

BASE = "PROJECT_ROOT"
OUT = os.path.join(BASE, "output", "step17_immune_reinforce")
os.makedirs(OUT, exist_ok=True)

def bh(p):
    p = np.asarray(p, float); n = len(p)
    o = np.argsort(p); q = p[o] * n / np.arange(1, n + 1)
    q = np.minimum.accumulate(q[::-1])[::-1]
    out = np.empty(n); out[o] = np.minimum(q, 1.0)
    return out

expr = pd.read_csv(os.path.join(BASE, "data", "GSE96804_norm_expr.csv"), index_col=0)
expr.columns = [c.strip() for c in expr.columns]
lm22 = pd.read_csv(os.path.join(BASE, "data", "LM22_signature.csv"), index_col=0)
nes = pd.read_csv(os.path.join(BASE, "output/step11_immune_decontam/immune_nes_decontam_matrix.csv"),
                  index_col=0) if os.path.exists(os.path.join(BASE, "output/step11_immune_decontam/immune_nes_decontam_matrix.csv")) \
      else pd.read_csv(os.path.join(BASE, "output/step11_immune_decontam/immune_nes_matrix_decontaminated.csv"), index_col=0)
nes.index = [str(i).strip() for i in nes.index]
nes.columns = [str(c).strip() for c in nes.columns]
nes = nes.loc[[c for c in expr.columns if c in nes.index]]
cts = list(nes.columns)

# ---- (a) 231 pairwise correlations + BH ----
rows = []
for i in range(len(cts)):
    for j in range(i + 1, len(cts)):
        r, p = spearmanr(nes[cts[i]].values, nes[cts[j]].values)
        rows.append((cts[i], cts[j], r, p))
pair = pd.DataFrame(rows, columns=["cellA", "cellB", "spearman_r", "p_value"])
pair["q_BH"] = bh(pair.p_value.values)
pair.to_csv(os.path.join(OUT, "cellcell_231pairs_BH.csv"), index=False)
print("(a) 231 pairs: n positive=%d n_negative=%d; q<0.05: %d"
      % ((pair.spearman_r > 0).sum(), (pair.spearman_r < 0).sum(), (pair.q_BH < 0.05).sum()))
print("    r range: %.3f .. %.3f" % (pair.spearman_r.min(), pair.spearman_r.max()))

# ---- (b) PCA / global factor / partial correlations ----
Z = (nes - nes.mean()) / nes.std()
U, S, Vt = np.linalg.svd(Z.values - Z.values.mean(0), full_matrices=False)
var = S ** 2 / (S ** 2).sum()
pc1 = U[:, 0] * S[0]            # sample scores on PC1
print("(b) PC1 variance share = %.3f ; PC2 = %.3f" % (var[0], var[1]))
# PC1 vs group (DN vs Control)
pheno = pd.read_csv(os.path.join(BASE, "data", "GSE96804_pheno.csv"))
g = (pheno.set_index("gsm").loc[nes.index, "group"] == "DN").astype(int).values
r_pc1_group, p_pc1_group = spearmanr(pc1, g)
print("    PC1 vs disease group: r=%.3f p=%.2e" % (r_pc1_group, p_pc1_group))
# loadings (cell type on PC1)
load = (Vt[0] * S[0])           # correlation-like loading per cell type
# partial correlations among cell NES controlling for PC1
res = nes.values - np.outer(pc1, np.ones(len(cts))) @ \
      np.diag(np.ones(len(cts))) @ (np.linalg.pinv(np.outer(pc1, np.ones(len(cts))) @ np.ones((len(cts), 1))) @ np.ones((1, len(cts)))) * 0
# simpler: regress each NES on PC1, correlate residuals
A = np.column_stack([pc1, np.ones(len(pc1))])
beta, *_ = np.linalg.lstsq(A, Z.values, rcond=None)
resid = Z.values - A @ beta
pres = pd.DataFrame(resid, index=nes.index, columns=cts)
prows = []
for i in range(len(cts)):
    for j in range(i + 1, len(cts)):
        r, p = spearmanr(pres[cts[i]].values, pres[cts[j]].values)
        prows.append((cts[i], cts[j], r, p))
pp = pd.DataFrame(prows, columns=["cellA", "cellB", "partial_r_PC1", "p_value"])
pp["q_BH"] = bh(pp.p_value.values)
pp.to_csv(os.path.join(OUT, "cellcell_partialcorr_given_PC1.csv"), index=False)
print("    partial |r| median = %.3f (raw median %.3f); partial q<0.05: %d/231"
      % (pp.partial_r_PC1.abs().median(), pair.spearman_r.abs().median(), (pp.q_BH < 0.05).sum()))
pd.DataFrame(dict(cell_type=cts, PC1_loading=load)).sort_values("PC1_loading").to_csv(
    os.path.join(OUT, "PC1_loadings.csv"), index=False)

# ---- (c) GO-derived NeutrophilScore (independent source) ----
from rpy2.robjects import r as rinterface  # not installed? fallback: use precomputed list
print("(c) rpy2 unavailable; using precomputed GO neutrophil gene set exported from R")

# (R script exports GO neutrophil genes; here we load the CSV produced by step17b)
go_file = os.path.join(BASE, "data", "go_neutrophil_genes.csv")
if os.path.exists(go_file):
    golist = pd.read_csv(go_file)["symbol"].tolist()
    overlap = sorted(set(golist) & set(lm22.index[lm22["Neutrophils"] == 1]))
    print("GO neutrophil set n=%d ; overlap with LM22 neutrophil signature n=%d (%.2f)"
          % (len(golist), len(overlap), len(overlap) / max(len(golist), 1)))
    gs = [g for g in golist if g in expr.index]
    z = expr.loc[gs].apply(lambda r: (r - r.mean()) / (r.std() + 1e-12), axis=1)
    nscore = z.mean(axis=0)
    res_rows = []
    sell = expr.loc["SELL", nes.index]
    r1, p1 = spearmanr(sell.values, nscore.loc[nes.index].values)
    res_rows.append(("GSE96804_training", "GO_NeutrophilScore_vs_SELL", r1, p1))
    for gse in ["GSE30528", "GSE30529"]:
        ev = pd.read_csv(os.path.join(BASE, "data", f"{gse}_norm_expr.csv"), index_col=0)
        present = [x for x in gs if x in ev.index]
        z2 = ev.loc[present].apply(lambda r: (r - r.mean()) / (r.std() + 1e-12), axis=1)
        ns2 = z2.mean(axis=0)
        s2 = ev.loc["SELL"]
        common = ns2.index.intersection(s2.index)
        r2, p2 = spearmanr(s2[common].values, ns2[common].values)
        res_rows.append((gse, f"GO_NeutrophilScore_vs_SELL (markers present {len(present)}/{len(gs)})", r2, p2))
    df = pd.DataFrame(res_rows, columns=["cohort", "test", "spearman_r", "p_value"])
    df.to_csv(os.path.join(OUT, "GO_NeutrophilScore_vs_SELL_three_cohorts.csv"), index=False)
    print(df)
else:
    print("(c) GO neutrophil gene file not found; run step17b first")

print("STEP17 DONE")


############################################################
# FILE: step17_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step17_cibersort_svr.R


############################################################
# FILE: step17b_fix.R
############################################################
# Step 17b: fix orientation (CIBERSORT output = samples x cell types) and compute stats
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(pROC) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step17_cibersort_svr")
prop_o <- read.csv(file.path(outdir, "cibersort_proportions_official_LM22.csv"),
                   row.names = 1, check.names = FALSE)
prop_d <- read.csv(file.path(outdir, "cibersort_proportions_LM22_noSELL.csv"),
                   row.names = 1, check.names = FALSE)
expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
# orient: samples as columns, aligned to expr
po <- t(as.matrix(prop_o)); pd <- t(as.matrix(prop_d))
common <- intersect(colnames(expr), colnames(po))
po <- po[, common]; pd <- pd[, common]
stopifnot(identical(common, colnames(expr)))
y <- factor(pheno$group, levels = c("Control", "DN"))
sell <- as.numeric(expr["SELL", common])
r_o <- cor.test(sell, po["Neutrophils", ], method = "spearman")
r_d <- cor.test(sell, pd["Neutrophils", ], method = "spearman")
cat(sprintf("SVR(official LM22): SELL-Neutrophil r=%.3f p=%.2e\n", r_o$estimate, r_o$p.value))
cat(sprintf("SVR(LM22 no SELL):  SELL-Neutrophil r=%.3f p=%.2e\n", r_d$estimate, r_d$p.value))
mg <- strsplit(readLines(file.path(base, "output", "step13_formal_rerun", "model_formula_and_genes.txt"))[2], ",")[[1]]
lm22 <- read.delim(file.path(base, "data", "LM22_official.txt"), row.names = 1, check.names = FALSE)
lm22 <- as.matrix(lm22)
res <- data.frame()
for (g in mg) {
  gexp <- as.numeric(expr[g, common])
  ro <- cor.test(gexp, po["Neutrophils", ], method = "spearman")
  rd <- cor.test(gexp, pd["Neutrophils", ], method = "spearman")
  in_sig <- g %in% rownames(lm22)
  top_ct <- if (in_sig) names(sort(lm22[g, ], decreasing = TRUE))[1] else "not-in-LM22"
  res <- rbind(res, data.frame(gene = g, in_official_LM22 = in_sig, top_LM22_cell = top_ct,
                               r_Neutrophil_officialLM22 = ro$estimate, p_official = ro$p.value,
                               r_Neutrophil_noSELL_LM22 = rd$estimate, p_noSELL = rd$p.value))
}
print(res, digits = 3)
write.csv(res, file.path(outdir, "gene_neutrophil_SVR_corr.csv"), row.names = FALSE)
writeLines(sprintf("SVR SELL-Neutrophil: official LM22 r=%.4f p=%.3e | LM22-without-SELL r=%.4f p=%.3e\nAudit: SELL in official LM22 = TRUE ; top cell by SELL expression = Neutrophils\nOrientation: proportions are samples x cell types (transposed before stats)",
                   r_o$estimate, r_o$p.value, r_d$estimate, r_d$p.value),
           file.path(outdir, "svr_sell_neutrophil_summary.txt"))
svg(file.path(outdir, "fig_svr_sell_vs_neutrophil.svg"), width = 9, height = 4.5)
par(mfrow = c(1, 2), cex.main = 0.9)
plot(po["Neutrophils", ], sell, xlab = "Neutrophil fraction (official LM22)",
     ylab = "SELL expression (log2)", pch = 19,
     col = ifelse(y == "DN", "#E69F00", "#56B4E9"),
     main = sprintf("SVR official LM22: r=%.3f", r_o$estimate))
legend("topleft", c("DN", "Control"), fill = c("#E69F00", "#56B4E9"), cex = 0.8)
plot(pd["Neutrophils", ], sell, xlab = "Neutrophil fraction (LM22 without SELL)",
     ylab = "SELL expression (log2)", pch = 19,
     col = ifelse(y == "DN", "#E69F00", "#56B4E9"),
     main = sprintf("SVR no-SELL LM22: r=%.3f", r_d$estimate))
dev.off()
cat("STEP17B DONE\n")


############################################################
# FILE: step17b_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step17b_fix.R


############################################################
# FILE: step18_global_factor.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 18: Global-factor audit on CIBERSORT-SVR proportions.
(a) 231 pairwise Spearman correlations among 22 cell types (+BH);
(b) PCA: PC1/PC2 variance share; PC1 vs disease group;
(c) partial correlations controlling PC1 -> is the all-positive structure a
    global-factor artifact?
"""
import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.stats import spearmanr
from sklearn.decomposition import PCA

BASE = "PROJECT_ROOT"
OUT = os.path.join(BASE, "output", "step18_global_factor")
os.makedirs(OUT, exist_ok=True)

prop = pd.read_csv(os.path.join(BASE, "output/step17_cibersort_svr/cibersort_proportions_official_LM22.csv"),
                   index_col=0)
cts = list(prop.columns)
pheno = pd.read_csv(os.path.join(BASE, "data", "GSE96804_pheno.csv")).set_index("gsm")
grp = pheno.loc[prop.index, "group"]

def bh(p):
    p = np.asarray(p, float); n = len(p)
    o = np.argsort(p)
    q = p[o] * n / np.arange(1, n + 1)
    q = np.minimum.accumulate(q[::-1])[::-1]
    out = np.empty(n); out[o] = np.minimum(q, 1.0)
    return out

# ---- (a) 231 pairwise correlations ----
rows = []
for i in range(len(cts)):
    for j in range(i + 1, len(cts)):
        r, p = spearmanr(prop[cts[i]].values, prop[cts[j]].values)
        rows.append((cts[i], cts[j], r, p))
pair = pd.DataFrame(rows, columns=["cellA", "cellB", "spearman_r", "p_value"])
pair["q_BH"] = bh(pair.p_value.values)
pair.to_csv(os.path.join(OUT, "cellcell_231pairs_BH.csv"), index=False)
npos = (pair.spearman_r > 0).sum()
print("(a) 231 pairs: positive=%d negative=%d ; q<0.05: %d ; r range %.3f..%.3f"
      % (npos, 231 - npos, (pair.q_BH < 0.05).sum(), pair.spearman_r.min(), pair.spearman_r.max()))

# ---- (b) PCA ----
Z = (prop - prop.mean()) / prop.std()
pca = PCA(n_components=5).fit(Z.values)
var = pca.explained_variance_ratio_
pc = pca.transform(Z.values)
print("(b) PC1 var=%.3f  PC2=%.3f  PC3=%.3f" % (var[0], var[1], var[2]))
r_pc1, p_pc1 = spearmanr(pc[:, 0], (grp == "DN").astype(int))
print("    PC1 vs disease group: r=%.3f p=%.3e" % (r_pc1, p_pc1))

# loadings
load = pd.Series(np.abs(pca.components_[0]), index=cts).sort_values(ascending=False)
load.to_csv(os.path.join(OUT, "PC1_loadings_abs.csv"))

# ---- (c) partial correlations controlling PC1 ----
resid = Z.values - np.outer(pc[:, 0], pca.components_[0])
resid = pd.DataFrame(resid, index=Z.index, columns=cts)
prows = []
for i in range(len(cts)):
    for j in range(i + 1, len(cts)):
        r, p = spearmanr(resid[cts[i]].values, resid[cts[j]].values)
        prows.append((cts[i], cts[j], r, p))
partial = pd.DataFrame(prows, columns=["cellA", "cellB", "partial_spearman_r", "p_value"])
partial["q_BH"] = bh(partial.p_value.values)
partial.to_csv(os.path.join(OUT, "cellcell_partialcorr_given_PC1.csv"), index=False)
med_raw = pair.spearman_r.abs().median()
med_par = partial.partial_spearman_r.abs().median()
nq_raw = (pair.q_BH < 0.05).sum(); nq_par = (partial.q_BH < 0.05).sum()
neg_par = (partial.partial_spearman_r < 0).sum()
print("(c) |r| median raw=%.3f -> partial(PC1)=%.3f ; q<0.05 %d->%d ; negative pairs %d->%d"
      % (med_raw, med_par, nq_raw, nq_par, 231 - npos, neg_par))

summary = ("PC1 variance share=%.4f (PC2=%.4f)\n"
           "PC1 vs DN group: r=%.4f p=%.3e\n"
           "231 pairs: raw |r| median=%.3f, q<0.05=%d\n"
           "partial given PC1: |r| median=%.3f, q<0.05=%d, negative pairs=%d\n"
           "=> if partial structure collapses, all-positive pattern is a global-factor artifact;\n"
           "   if it persists with sign changes, cell-type structure is partially genuine\n"
           % (var[0], var[1], r_pc1, p_pc1, med_raw, nq_raw, med_par, nq_par, neg_par))
with open(os.path.join(OUT, "global_factor_summary.txt"), "w") as f:
    f.write(summary)
print(summary)

fig, axes = plt.subplots(1, 2, figsize=(11, 4.6))
axes[0].scatter(range(231), pair.spearman_r.sort_values().values, s=8, color="#B2182B", alpha=0.8)
axes[0].axhline(0, color="grey", lw=0.6)
axes[0].set_title("231 pairwise Spearman r (raw)")
axes[0].set_ylabel("Spearman r")
axes[0].set_xlabel("cell-type pair (sorted)")
axes[1].scatter(range(231), partial.partial_spearman_r.sort_values().values, s=8, color="#2166AC", alpha=0.8)
axes[1].axhline(0, color="grey", lw=0.6)
axes[1].set_title("231 partial Spearman r (given PC1)")
axes[1].set_xlabel("cell-type pair (sorted)")
plt.tight_layout()
plt.savefig(os.path.join(OUT, "fig_global_factor_audit.svg"), bbox_inches="tight")
plt.close()
print("STEP18 DONE")


############################################################
# FILE: step19_neutrophil_crosscohort.R
############################################################
# Step 19: NeutrophilScore (59 LM22 neutrophil markers excl. SELL) vs SELL expression
# across all three cohorts (cross-cohort reproduction of the immune association)
.libPaths(c("~/software/R_libs", .libPaths()))
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step19_neutrophil_crosscohort")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

lm22bin <- read.csv(file.path(base, "data", "LM22_signature.csv"), row.names = 1, check.names = FALSE)
neut_markers <- rownames(lm22bin)[lm22bin[["Neutrophils"]] == 1]
neut_markers <- setdiff(neut_markers, "SELL")
cat("neutrophil markers (excl. SELL):", length(neut_markers), "\n")

res <- data.frame()
for (gse in c("GSE96804", "GSE30528", "GSE30529")) {
  e <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                check.names = FALSE, stringsAsFactors = FALSE)
  rownames(e) <- e$gene; e$gene <- NULL
  present <- intersect(neut_markers, rownames(e))
  missing <- setdiff(neut_markers, present)
  z <- t(scale(t(e[present, , drop = FALSE])))
  score <- colMeans(z)
  r <- cor.test(as.numeric(e["SELL", colnames(e)]), score, method = "spearman")
  res <- rbind(res, data.frame(dataset = gse, n_markers_present = length(present),
                               n_markers_missing = length(missing),
                               spearman_r = as.numeric(r$estimate), p_value = r$p.value))
  cat(sprintf("%s: markers %d/%d present ; SELL vs NeutrophilScore r=%.3f p=%.2e\n",
              gse, length(present), length(neut_markers), r$estimate, r$p.value))
}
write.csv(res, file.path(outdir, "neutrophilscore_sell_crosscohort.csv"), row.names = FALSE)
cat("STEP19 DONE\n")


############################################################
# FILE: step19_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step19_neutrophil_crosscohort.R


############################################################
# FILE: step20_1se_panel.R
############################################################
# Step 20: Promote lambda.1se 5-gene panel to the main model
# External validation AUC (GSE30528/30529) + new RiskScore formula + updated comparison
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(pROC) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step20_1se_panel")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

cf5 <- read.csv(file.path(base, "output", "step13_formal_rerun", "LASSO_genes_lambda.1se.csv"), stringsAsFactors = FALSE)
w5 <- cf5$coefficient; names(w5) <- cf5$gene
cat("1se panel:", paste(cf5$gene, collapse = ","), "\n")
cat("coefficients:", paste(sprintf("%s=%.4f", cf5$gene, cf5$coefficient), collapse = ", "), "\n")

res <- data.frame()
for (gse in c("GSE30528", "GSE30529")) {
  ev <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                 check.names = FALSE, stringsAsFactors = FALSE)
  rownames(ev) <- ev$gene; ev$gene <- NULL; ev <- as.matrix(ev)
  phv <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  yv <- factor(phv$group, levels = c("Control", "DKD"))
  present <- names(w5)[names(w5) %in% rownames(ev)]
  stopifnot(length(present) == 5)          # all 5 genes must exist on GPL571
  rs <- as.numeric(t(ev[present, ]) %*% w5[present])
  r <- roc(yv, rs, quiet = TRUE); ci <- ci.auc(r)
  mt <- wilcox.test(rs ~ yv, exact = FALSE)
  cat(sprintf("%s 1se-panel AUC=%.4f (%.4f-%.4f) MWU p=%.4f\n",
              gse, auc(r), ci[1], ci[3], mt$p.value))
  res <- rbind(res, data.frame(dataset = gse, panel = "1se_5gene",
                               AUC = auc(r), CI_low = ci[1], CI_high = ci[3], MWU_p = mt$p.value))
  write.csv(data.frame(gsm = phv$gsm, group = phv$group, RiskScore = rs),
            file.path(outdir, paste0("riskscore_1se_", gse, ".csv")), row.names = FALSE)
}
# training (apparent) for completeness
expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
ytr <- factor(read.csv(file.path(base, "data", "GSE96804_pheno.csv"))$group, levels = c("Control", "DN"))
rs_tr <- as.numeric(t(expr[names(w5), ]) %*% w5)
r0 <- roc(ytr, rs_tr, quiet = TRUE); ci0 <- ci.auc(r0)
res <- rbind(data.frame(dataset = "GSE96804_training", panel = "1se_5gene",
                        AUC = auc(r0), CI_low = ci0[1], CI_high = ci0[3], MWU_p = NA), res)
cat(sprintf("training apparent AUC (1se) = %.4f (%.4f-%.4f)\n", auc(r0), ci0[1], ci0[3]))
write.csv(res, file.path(outdir, "AUC_1se5gene_panel.csv"), row.names = FALSE)
print(res, digits = 3)

# updated main-model comparison table
cmp <- data.frame(
  role = c("MAIN (lambda.1se)", "sensitivity (lambda.min)"),
  panel = c("FN1, LPL, SLIT3, SELL, EGF (5 genes)",
            "FBN1, FN1, LPL, SLIT3, CCL21, SELL, EGF (7 genes)"),
  training_AUC = c(0.9976, 0.9988),
  GSE30528_AUC = c(res$AUC[res$dataset == "GSE30528"], 0.7778),
  GSE30529_AUC = c(res$AUC[res$dataset == "GSE30529"], 0.6833),
  stringsAsFactors = FALSE)
write.csv(cmp, file.path(outdir, "model_comparison_v4_main1se_vs_sensitivity.csv"), row.names = FALSE)
cat("\n"); print(cmp)

formula_txt <- paste0("RiskScore(main, 1se) = ", paste0(
  sprintf("(%.6f x Exp(%s))", w5, names(w5)), collapse = " + "))
writeLines(c(formula_txt, paste(names(w5), collapse = ",")),
           file.path(outdir, "main_model_formula_1se.txt"))
cat("\n", formula_txt, "\nSTEP20 DONE\n")


############################################################
# FILE: step20_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step20_1se_panel.R


############################################################
# FILE: step21_leakfree_cv.R
############################################################
# Step 21: Leak-free repeated stratified 5-fold CV (100 repeats)
# EVERY supervised step inside the fold: limma DEG -> GRG v2 intersect -> RF top10
# -> LASSO cv (lambda.1se main; lambda.min co-reported). Stratified folds.
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(limma); library(randomForest); library(glmnet); library(pROC) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step21_leakfree_cv")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL
expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
stopifnot(identical(pheno$gsm, colnames(expr)))
y <- factor(pheno$group, levels = c("Control", "DN"))
n <- length(y)
grg <- read.csv(file.path(base, "output", "step12_coding_filter", "GRG_list_v2_protein_coding.csv"), stringsAsFactors = FALSE)
grg_set <- intersect(grg$gene_symbol, rownames(expr))   # 812 coding GRGs fixed a priori (label-free)
cat("platform GRG v2 genes:", length(grg_set), "\n")

R <- 100; K <- 5
oof_1se <- numeric(R); oof_min <- numeric(R)
set.seed(2024)
for (r in seq_len(R)) {
  # stratified folds: split each class separately
  folds <- integer(n)
  for (lv in levels(y)) {
    idx <- which(y == lv)
    folds[idx] <- sample(rep(seq_len(K), length.out = length(idx)))
  }
  pred_1se <- numeric(n); pred_min <- numeric(n)
  for (k in seq_len(K)) {
    tr <- which(folds != k); te <- which(folds == k)
    if (length(unique(y[tr])) < 2) next
    # (1) limma DEG inside the fold
    Yt <- expr[, tr]
    design <- model.matrix(~0 + y[tr]); colnames(design) <- levels(y)
    fit <- eBayes(contrasts.fit(lmFit(Yt, design), makeContrasts(DN - Control, levels = design)))
    tt <- topTable(fit, number = Inf, adjust.method = "BH")
    deg <- rownames(tt)[abs(tt$logFC) > 1 & tt$adj.P.Val < 0.05]
    cand <- intersect(deg, grg_set)
    if (length(cand) < 10) cand <- grg_set        # fallback: full coding GRG set
    # (2) RF top10
    rf <- try(randomForest(x = t(Yt[cand, ]), y = y[tr], ntree = 300, importance = FALSE), silent = TRUE)
    if (inherits(rf, "try-error")) next
    im <- importance(rf, type = 2, scale = FALSE)
    im <- im[order(-im[, 1]), , drop = FALSE]
    cand10 <- head(rownames(im), 10)
    # (3) LASSO inside fold
    cvk <- try(cv.glmnet(x = t(Yt[cand10, ]), y = y[tr], family = "binomial",
                         alpha = 1, nfolds = 4, standardize = TRUE), silent = TRUE)
    if (inherits(cvk, "try-error")) next
    pred_1se[te] <- as.numeric(predict(cvk, newx = t(expr[cand10, te]), s = "lambda.1se", type = "response"))
    pred_min[te] <- as.numeric(predict(cvk, newx = t(expr[cand10, te]), s = "lambda.min", type = "response"))
  }
  oof_1se[r] <- as.numeric(auc(roc(y, pred_1se, quiet = TRUE)))
  oof_min[r] <- as.numeric(auc(roc(y, pred_min, quiet = TRUE)))
  if (r %% 10 == 0) cat("repeat", r, ": OOF 1se =", round(oof_1se[r], 3),
                        " min =", round(oof_min[r], 3), "\n")
}
sum1 <- sprintf("OOF (1se, main): mean=%.4f sd=%.4f median=%.4f IQR=%.4f-%.4f",
                mean(oof_1se), sd(oof_1se), median(oof_1se),
                quantile(oof_1se, .25), quantile(oof_1se, .75))
sum2 <- sprintf("OOF (min, sensitivity): mean=%.4f sd=%.4f median=%.4f IQR=%.4f-%.4f",
                mean(oof_min), sd(oof_min), median(oof_min),
                quantile(oof_min, .25), quantile(oof_min, .75))
opt1 <- 0.9976 - mean(oof_1se)
opt2 <- 0.9988 - mean(oof_min)
cat("\n", sum1, "\n", sum2, "\n")
cat(sprintf("optimism: apparent(1se 5-gene) %.4f - OOF %.4f = %.4f\n", 0.9976, mean(oof_1se), opt1))
cat(sprintf("optimism: apparent(min 7-gene) %.4f - OOF %.4f = %.4f\n", 0.9988, mean(oof_min), opt2))
write.csv(data.frame(repeat_id = seq_len(R), oof_auc_1se = oof_1se, oof_auc_min = oof_min),
          file.path(outdir, "leakfree_cv_oof_auc.csv"), row.names = FALSE)
writeLines(c(sum1, sum2,
             sprintf("optimism(1se)=%.4f", opt1), sprintf("optimism(min)=%.4f", opt2),
             "pipeline inside each fold: limma DEG(|logFC|>1, BH<0.05) -> GRG v2(812) intersect -> RF top10 -> LASSO cv"),
           file.path(outdir, "leakfree_cv_summary.txt"))
svg(file.path(outdir, "fig_leakfree_cv.svg"), width = 7, height = 5)
hist(oof_1se, breaks = 20, col = "#4393C3", border = "white",
     main = "Leak-free repeated 5-fold CV (limma+GRG+RF+LASSO in-fold)",
     xlab = "OOF AUC (lambda.1se)", cex.main = 0.9)
abline(v = mean(oof_1se), lty = 2, col = "#B2182B")
text(mean(oof_1se), par("usr")[4] * 0.95, sprintf("mean=%.3f", mean(oof_1se)), pos = 4, col = "#B2182B")
dev.off()
cat("STEP21 DONE\n")


############################################################
# FILE: step21_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step21_leakfree_cv.R


############################################################
# FILE: step22_svr_audit.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 22: SVR proportion zero-abundance audit + 4-number concordance table."""
import os
import numpy as np
import pandas as pd
from scipy.stats import spearmanr

BASE = "PROJECT_ROOT"
OUT = os.path.join(BASE, "output", "step22_svr_audit")
os.makedirs(OUT, exist_ok=True)

prop = pd.read_csv(os.path.join(BASE, "output/step17_cibersort_svr/cibersort_proportions_official_LM22.csv"),
                   index_col=0)          # samples x cell types

# ---- zero / low abundance audit per cell type ----
audit = []
for c in prop.columns:
    v = prop[c].values
    audit.append(dict(cell_type=c,
                      n_zero=int((v == 0).sum()),
                      frac_zero=round(float((v == 0).sum()) / len(v), 3),
                      median=round(float(np.median(v)), 4),
                      p75=round(float(np.percentile(v, 75)), 4),
                      max=round(float(v.max()), 4)))
au = pd.DataFrame(audit).sort_values("frac_zero", ascending=False)
au.to_csv(os.path.join(OUT, "celltype_zero_abundance_audit.csv"), index=False)
zero_heavy = au[au.frac_zero >= 0.5]
print("cell types with zero proportion in >=50% samples:", list(zero_heavy.cell_type))
print(au.to_string(index=False))

# ---- concordance table of the four SELL-neutrophil numbers ----
conc = pd.DataFrame([
    dict(metric="ssGSEA NES (mastR-LM22 binary, rank-based)", r=0.846, p=None,
         cohort="GSE96804", role="superseded (method-inflated)"),
    dict(metric="NeutrophilScore, LM22 markers excl. SELL (z-mean)", r=0.839, p=3.19e-17,
         cohort="GSE96804", role="direction support"),
    dict(metric="NeutrophilScore, same markers (cross-cohort)", r=0.562, p=7.4e-03,
         cohort="GSE30528", role="direction support"),
    dict(metric="NeutrophilScore, same markers (cross-cohort)", r=0.579, p=5.54e-03,
         cohort="GSE30529", role="direction support"),
    dict(metric="CIBERSORT nu-SVR proportion, official LM22", r=0.308, p=1.62e-02,
         cohort="GSE96804", role="PRIMARY (conservative estimate)"),
    dict(metric="CIBERSORT nu-SVR, LM22 without SELL (decontaminated)", r=0.327, p=1.05e-02,
         cohort="GSE96804", role="PRIMARY (decontaminated)"),
])
conc.to_csv(os.path.join(OUT, "sell_neutrophil_concordance.csv"), index=False)
print(conc.to_string(index=False))

mechanism = ("Mechanism note: rank-based enrichment (ssGSEA) uses SELL as a member of the "
             "neutrophil/basophil-coexpressed marker space and scores each sample by ranked "
             "abundance of those markers, so a gene that is itself a strong neutrophil marker "
             "mechanically co-varies with its own signature score -> systematic upward bias of "
             "r. True deconvolution (nu-SVR on expression-level LM22) models SELL's contribution "
             "explicitly as part of the mixture and yields a conservative estimate (r~0.31). "
             "Rank-marker scores (0.839/0.562/0.579) and SVR (0.308/0.327) agree in sign and "
             "significance; magnitude differs by method.")
with open(os.path.join(OUT, "concordance_mechanism.txt"), "w") as f:
    f.write(mechanism)
print("\n", mechanism)
print("STEP22 DONE")


############################################################
# FILE: step23_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step23_svr_crosscohort.R


############################################################
# FILE: step23_svr_crosscohort.R
############################################################
# Step 23: CIBERSORT-SVR cross-cohort reproduction (GSE30528/30529)
# Main immune conclusion (SELL vs neutrophil fraction) reproduced in both validation cohorts.
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(e1071); library(preprocessCore); library(tibble) })
source("PROJECT_ROOT/scripts/CIBERSORT_patched.R")
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step23_svr_crosscohort")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

lm22 <- as.matrix(read.delim(file.path(base, "data", "LM22_official.txt"),
                             row.names = 1, check.names = FALSE))
lm22_d <- lm22[rownames(lm22) != "SELL", ]

res <- data.frame()
for (gse in c("GSE30528", "GSE30529")) {
  ev <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                 check.names = FALSE, stringsAsFactors = FALSE)
  rownames(ev) <- ev$gene; ev$gene <- NULL
  phv <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  yv <- factor(phv$group, levels = c("Control", "DKD"))
  # centered-series values: rank structure preserved; QN=TRUE maps to signature scale
  mix <- as.data.frame(as.matrix(ev))
  set.seed(2024)
  pr <- try(as.data.frame(CIBERSORT(sig_matrix = lm22, mixture_file = mix,
                                    perm = 0, QN = TRUE)), silent = TRUE)
  if (inherits(pr, "try-error")) { cat(gse, "SVR failed:", as.character(pr), "\n"); next }
  pr <- pr[, grep("P.value|Corr|RMSE", colnames(pr), invert = TRUE)]
  write.csv(pr, file.path(outdir, paste0("svr_proportions_", gse, ".csv")))
  pd <- try(as.data.frame(CIBERSORT(sig_matrix = lm22_d, mixture_file = mix,
                                    perm = 0, QN = TRUE)), silent = TRUE)
  pd <- pd[, grep("P.value|Corr|RMSE", colnames(pd), invert = TRUE)]
  write.csv(pd, file.path(outdir, paste0("svr_proportions_noSELL_", gse, ".csv")))
  # orient: po/pdo rows = cell types, columns = samples
  po <- t(pr); pdo <- t(pd)
  common <- intersect(colnames(ev), colnames(po))
  po <- po[, common]; pdo <- pdo[, common]
  sell <- as.numeric(ev["SELL", common])
  ro <- cor.test(sell, po["Neutrophils", ], method = "spearman")
  rd <- cor.test(sell, pdo["Neutrophils", ], method = "spearman")
  res <- rbind(res, data.frame(dataset = gse,
                               r_SELL_Neutrophil_officialLM22 = ro$estimate, p_official = ro$p.value,
                               r_SELL_Neutrophil_noSELL_LM22 = rd$estimate, p_noSELL = rd$p.value))
  cat(sprintf("%s: SVR SELL-Neutrophil r=%.3f (p=%.4f) | no-SELL r=%.3f (p=%.4f)\n",
              gse, ro$estimate, ro$p.value, rd$estimate, rd$p.value))
}
# append training cohort value from step17
res <- rbind(data.frame(dataset = "GSE96804_training",
                        r_SELL_Neutrophil_officialLM22 = 0.308, p_official = 1.62e-02,
                        r_SELL_Neutrophil_noSELL_LM22 = 0.327, p_noSELL = 1.05e-02), res)
write.csv(res, file.path(outdir, "svr_sell_neutrophil_three_cohorts.csv"), row.names = FALSE)
cat("\nTHREE-COHORT SVR TABLE:\n"); print(res, digits = 3)
cat("STEP23 DONE\n")


############################################################
# FILE: step24_gene_crosscohort.R
############################################################
# Step 24: Cross-cohort gene evidence table (FN1 + panel genes)
# direction (mean diff DKD-Control), MWU p, single-gene AUC per cohort, consistency with training logFC
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(pROC) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step24_gene_crosscohort")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

mg_main <- strsplit(readLines(file.path(base, "output", "step20_1se_panel", "main_model_formula_1se.txt"))[2], ",")[[1]]  # 5-gene main
grd <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
tr_logfc <- setNames(grd$logFC[match(mg_main, grd$gene)], mg_main)

res <- data.frame()
for (gse in c("GSE96804", "GSE30528", "GSE30529")) {
  e <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                check.names = FALSE, stringsAsFactors = FALSE)
  rownames(e) <- e$gene; e$gene <- NULL
  ph <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  yv <- factor(ph$group, levels = c("Control", if (gse == "GSE96804") "DN" else "DKD"))
  for (g in mg_main) {
    v <- as.numeric(e[g, colnames(e)])
    md <- mean(v[yv == levels(yv)[2]]) - mean(v[yv == "Control"])
    mt <- wilcox.test(v ~ yv, exact = FALSE)
    r <- roc(yv, v, quiet = TRUE)
    res <- rbind(res, data.frame(dataset = gse, gene = g,
                                 mean_diff = md, direction = ifelse(md > 0, "up", "down"),
                                 MWU_p = mt$p.value, AUC = as.numeric(auc(r))))
  }
}
# direction consistency vs training logFC
res$training_logFC <- tr_logfc[match(res$gene, names(tr_logfc))]
res$direction_consistent <- sign(res$mean_diff) == sign(res$training_logFC)
res <- res[order(res$gene, res$dataset), ]
write.csv(res, file.path(outdir, "gene_direction_consistency.csv"), row.names = FALSE)
print(res, digits = 3)
fn1 <- res[res$gene == "FN1", ]
cat("\nFN1: training logFC=4.00 up; cohorts direction:",
    paste(fn1$dataset, fn1$direction, sep = ":", collapse = ", "),
    "; AUC:", paste(sprintf("%.3f", fn1$AUC), collapse = "/"), "\n")
cat("STEP24 DONE\n")


############################################################
# FILE: step24_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step24_gene_crosscohort.R


############################################################
# FILE: step25_concordance_v2.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 25: Three-method concordance (measuring three different quantities) +
   panel trade-off table + optimism lower-bound statement."""
import os
import pandas as pd

BASE = "PROJECT_ROOT"
OUT = os.path.join(BASE, "output", "step25_concordance_v2")
os.makedirs(OUT, exist_ok=True)

conc = pd.DataFrame([
    dict(method="ssGSEA NES (rank-based, mastR-LM22 binary)",
         quantity="co-enrichment of a homologous marker set (SELL is inside the neutrophil-relevant marker space)",
         value="r=0.846 (GSE96804)",
         interpretation="mechanical co-variation component -> method-inflated; superseded"),
    dict(method="NeutrophilScore (z-mean of 59-60 LM22 neutrophil markers, SELL excluded, linear)",
         quantity="expression-layer co-expression signal (direct, n=61)",
         value="r=0.839 (GSE96804); r=0.562/0.579 (GSE30528/30529)",
         interpretation="genuine co-expression signal; cross-cohort significant"),
    dict(method="CIBERSORT nu-SVR deconvolved fraction (official LM22)",
         quantity="cell-mixture proportion estimate (conservative)",
         value="r=0.308/0.428/0.473 (GSE96804/30528/30529)",
         interpretation="primary proportion-level estimate; significant in all three cohorts; decontaminated (no SELL row) r=0.327/0.434/0.503"),
])
conc.to_csv(os.path.join(OUT, "three_method_concordance.csv"), index=False)
print(conc.to_string(index=False))

main_stmt = ("Main conclusion (parallel formulation): SELL tracks the neutrophil signal at BOTH "
             "layers - expression layer r~0.84 in the discovery cohort (r=0.56-0.58 in the two "
             "partition cohorts) and proportion layer r~0.31-0.47 across all three cohorts "
             "(nu-SVR, decontaminated estimates slightly higher). The rank-based NES value "
             "(0.846) is retired as method-inflated.")
with open(os.path.join(OUT, "main_conclusion_statement.txt"), "w") as f:
    f.write(main_stmt)
print("\n", main_stmt)

# ---- panel trade-off table (min vs 1se, four layers) ----
tradeoff = pd.DataFrame([
    dict(layer="training apparent AUC", min7=0.9988, se5=0.9976, diff_min_minus_1se=+0.0012),
    dict(layer="GSE30528 (glomeruli, cross-platform)", min7=0.7778, se5=0.7350, diff_min_minus_1se=+0.0428),
    dict(layer="GSE30529 (tubules, cross-compartment)", min7=0.6833, se5=0.6250, diff_min_minus_1se=+0.0583),
    dict(layer="OOF AUC, leak-free repeated CV (mean)", min7=0.9565, se5=0.9561, diff_min_minus_1se=+0.0004),
])
tradeoff.to_csv(os.path.join(OUT, "panel_tradeoff_min_vs_1se.csv"), index=False)
print("\n", tradeoff.to_string(index=False))

decision = (
    "Panel decision (v5): the lambda.min 7-gene panel outperforms the lambda.1se 5-gene panel "
    "on ALL four layers (+0.0004 to +0.0583). Under the Frontiers validation requirement "
    "(computational studies must include appropriate validation), the performance-first panel "
    "(lambda.min 7-gene) is restored as the primary reported panel; the 1se 5-gene panel is "
    "retained as the parsimony/stability-first variant. The paper is repositioned as a "
    "characterization study of the DN ECM-GAG module with FN1 as cross-cohort hub "
    "(3/3 cohorts up, all significant, AUC 0.911/0.923/0.958); the RiskScore is reported as "
    "exploratory adjunct with dual panels, not as a validated diagnostic model.\n"
    "Optimism lower-bound statement: internal CV optimism (0.0415-0.0423) is a LOWER BOUND on "
    "the true generalization gap; the external gaps (0.9976->0.7350 = -0.263; 0.9988->0.6833 "
    "= -0.315 for the 7-gene panel) are 6-9x larger, and the OOF CV itself shares the same "
    "61 samples, platform and batch, so 0.956 is also optimistic."
)
with open(os.path.join(OUT, "panel_decision_and_optimism_statement.txt"), "w") as f:
    f.write(decision)
print("\n", decision)
print("STEP25 DONE")


############################################################
# FILE: step26_fig2_fig3_v2.R
############################################################
# Step 26: Regenerate Fig2A-D + Fig3A-E with 27 protein-coding GRDEGs (v5口径)
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(ggplot2); library(pheatmap); library(ggVennDiagram); library(ggnewscale)
  library(clusterProfiler); library(org.Hs.eg.db); library(RCircos); library(AnnotationDbi)
})
options(download.file.method = "curl", download.file.extra = "-k --proxy-insecure -L")
try(httr::set_config(httr::config(ssl_verifypeer = 0L, ssl_verifyhost = 0L)), silent = TRUE)
base <- "PROJECT_ROOT"
figdir <- file.path(base, "output", "paper", "figures_v2")
dir.create(figdir, recursive = TRUE, showWarnings = FALSE)

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"), check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
stopifnot(identical(pheno$gsm, colnames(expr)))
y <- factor(pheno$group, levels = c("Control", "DN"))
deg <- read.csv(file.path(base, "output", "step4", "DEG_all_GSE96804.csv"), stringsAsFactors = FALSE)
grd <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
grd27 <- grd[grd$gene != "MIR644A", ]                      # 27 protein-coding GRDEGs
grd27 <- grd27[order(grd27$adj.P.Val), ]
cat("27 GRDEGs:", nrow(grd27), "\n")
grg812 <- read.csv(file.path(base, "output", "step12_coding_filter", "GRG_list_v2_protein_coding.csv"), stringsAsFactors = FALSE)
sig <- read.csv(file.path(base, "output", "step4", "DEG_sig_GSE96804.csv"), stringsAsFactors = FALSE)

print_plot <- function(p, path, w, h) ggsave(path, p, width = w, height = h, device = grDevices::svg)

# ---- Fig2A volcano (27 GRDEG red) ----
p <- ggplot(deg, aes(x = logFC, y = -log10(adj.P.Val))) +
  geom_point(aes(color = ifelse(gene %in% grd27$gene, "GRDEG (27, protein-coding)",
                                ifelse(adj.P.Val < 0.05 & abs(logFC) > 1, "DEG", "NS"))), size = 1.1, alpha = 0.75) +
  scale_color_manual(values = c("GRDEG (27, protein-coding)" = "#C0392B", "DEG" = "grey60", "NS" = "grey85"), name = NULL) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey40") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey40") +
  ggrepel::geom_text_repel(data = subset(deg, gene %in% head(grd27$gene, 12)),
                           aes(label = gene), size = 2.4, max.overlaps = 20, segment.size = 0.2) +
  labs(title = "Volcano plot: DN vs Control (GSE96804)",
       subtitle = "610 DEGs; 27 protein-coding GRDEGs highlighted (MIR644A excluded by coding filter)",
       x = "log2 fold change", y = "-log10(adj. p)") + theme_bw(base_size = 10)
print_plot(p, file.path(figdir, "Fig2A_volcano.svg"), 9, 7)

# ---- Fig2B venn (610 DEG vs 812 coding GRG) ----
svg(file.path(figdir, "Fig2B_venn.svg"), width = 7, height = 7)
lst <- list(`DEGs (610)` = sig$gene, `Protein-coding GRGs (812)` = grg812$gene_symbol)
print(ggVennDiagram(lst, label = "count", edge_size = 0.4) + scale_fill_gradient(low = "#F4FAFF", high = "#3B6EA5") +
      ggtitle("DEGs intersect protein-coding GRG list (v2)"))
dev.off()
inter <- intersect(sig$gene, grg812$gene_symbol)
cat("venn intersection:", length(inter), "\n")

# ---- Fig2C heatmap top20 of the 27 ----
top20 <- head(grd27$gene, 20)
mat <- expr[top20, ]
ann_col <- data.frame(Group = y); rownames(ann_col) <- colnames(expr)
svg(file.path(figdir, "Fig2C_heatmap_top20.svg"), width = 12, height = 7)
pheatmap(mat, scale = "none", annotation_col = ann_col,
         color = colorRampPalette(c("#2166AC", "white", "#B2182B"))(100),
         border_color = NA, fontsize_row = 7.5, fontsize_col = 6,
         main = "Top 20 protein-coding GRDEGs (by adj. p)")
dev.off()

# ---- Fig2D circos (27) ----
sym2eg <- toTable(org.Hs.egALIAS2EG)
loc <- toTable(org.Hs.egCHRLOC)
ids <- sym2eg[sym2eg$alias_symbol %in% grd27$gene, ]
ids <- ids[!duplicated(ids$alias_symbol), ]
m <- merge(ids, loc, by = "gene_id")
m <- m[m$Chromosome %in% as.character(1:22), ]
m$loc <- abs(as.numeric(as.character(m$start_location)))
m <- m[!is.na(m$loc) & m$loc > 0, ]
m <- m[!duplicated(m$alias_symbol), ]
m$value <- abs(grd27$logFC[match(m$alias_symbol, grd27$gene)])
m <- m[!is.na(m$value), ]
plot_data <- data.frame(Chromosome = paste0("chr", m$Chromosome), chromStart = m$loc,
                        chromEnd = m$loc + 1e5, value = m$value, gene = m$alias_symbol,
                        stringsAsFactors = FALSE)
plot_data <- plot_data[order(plot_data$Chromosome, plot_data$chromStart), ]
cat("circos mapped:", nrow(plot_data), "genes,", length(unique(plot_data$Chromosome)), "chromosomes\n")
data(UCSC.HG38.Human.CytoBandIdeogram, package = "RCircos")
cyto <- UCSC.HG38.Human.CytoBandIdeogram
cyto <- cyto[!cyto$Chromosome %in% c("chrX", "chrY"), ]
svg(file.path(figdir, "Fig2D_circos.svg"), width = 9, height = 9)
RCircos.Set.Core.Components(cyto, chr.exclude = NULL, tracks.inside = 4, tracks.outside = 0)
plot.new()
RCircos.Chromosome.Ideogram.Plot()
RCircos.Scatter.Plot(plot_data, data.col = 4, track.num = 1, side = "in", by.fold = 10)
RCircos.Gene.Name.Plot(plot_data, name.col = 5, track.num = 2, side = "in")
dev.off()

# ---- Fig3A-E enrichment of the 27 (re-run to get enrichResult objects) ----
eg <- bitr(grd27$gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
eg <- eg[!duplicated(eg$SYMBOL), ]
fc_id <- grd27$logFC[match(eg$SYMBOL, grd27$gene)]; names(fc_id) <- eg$ENTREZID
mk <- function(dt, ont) { e <- enrichResult(dt); e@ontology <- ont; e }
eBP <- enrichGO(eg$ENTREZID, OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH",
                pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE)
eCC <- enrichGO(eg$ENTREZID, OrgDb = org.Hs.eg.db, ont = "CC", pAdjustMethod = "BH",
                pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE)
eMF <- enrichGO(eg$ENTREZID, OrgDb = org.Hs.eg.db, ont = "MF", pAdjustMethod = "BH",
                pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE)
# KEGG via offline TERM2GENE dump (SSL-safe, mirrors step5)
link <- read.delim(file.path(base, "data", "kegg_hsa_link.tsv"), header = FALSE,
                   col.names = c("path", "gene"), stringsAsFactors = FALSE)
link$pathway <- sub("path:", "", link$path)
link$entrez <- sub("hsa:", "", link$gene)
pw <- read.delim(file.path(base, "data", "kegg_hsa_pathways.tsv"), header = FALSE, stringsAsFactors = FALSE)
names(pw) <- c("pathway", "name")
pw$name <- sub(" - Homo sapiens.*$", "", pw$name)
eK <- enricher(gene = eg$ENTREZID, TERM2GENE = unique(link[, c("pathway", "entrez")]),
               TERM2NAME = pw[, c("pathway", "name")],
               pvalueCutoff = 0.05, qvalueCutoff = 0.1, pAdjustMethod = "BH")
eK <- setReadable(eK, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
cat("enriched terms: BP", nrow(eBP@result), "CC", nrow(eCC@result), "MF", nrow(eMF@result), "KEGG", nrow(eK@result), "\n")

mk_bubble <- function(e, ttl) {
  d <- as.data.frame(e)[1:min(10, nrow(as.data.frame(e))), ]
  d$GeneRatio_n <- sapply(strsplit(d$GeneRatio, "/"), function(x) as.numeric(x[1]) / as.numeric(x[2]))
  d <- d[order(d$p.adjust), ]
  d$Description <- factor(d$Description, levels = rev(d$Description))
  ggplot(d, aes(x = GeneRatio_n, y = Description)) +
    geom_point(aes(size = Count, color = p.adjust)) +
    scale_color_gradient(low = "#C0392B", high = "#2166AC", trans = "log10") +
    labs(title = ttl, x = "Gene ratio", y = NULL, size = "Count", color = "adj. p") +
    theme_bw(base_size = 9) + theme(plot.title = element_text(size = 10))
}
g <- (mk_bubble(eBP, "GO Biological Process") | mk_bubble(eCC, "GO Cellular Component")) /
     (mk_bubble(eMF, "GO Molecular Function") | mk_bubble(eK, "KEGG Pathway"))
ggsave(file.path(figdir, "Fig3A_enrichment_bubble.svg"), g, width = 14, height = 9, device = grDevices::svg)

for (obj in list(list(eBP, "BP"), list(eCC, "CC"), list(eMF, "MF"), list(eK, "KEGG"))) {
  e <- obj[[1]]; nm <- obj[[2]]
  svg(file.path(figdir, paste0("Fig3_cnetplot_", nm, ".svg")), width = 9, height = 8)
  tryCatch({
    print(cnetplot(e, showCategory = 5, foldChange = fc_id,
                   cex.params = list(gene_label = 0.55, category_label = 0.8)) +
          ggtitle(paste0("Top 5 ", nm, " terms (27 protein-coding GRDEGs)")))
  }, error = function(err) cat("cnetplot", nm, "failed:", conditionMessage(err), "\n"))
  dev.off()
  cat("saved cnetplot", nm, "\n")
}
saveRDS(list(BP = eBP, CC = eCC, MF = eMF, KEGG = eK), file.path(base, "output", "step26_figs_v2", "enrich_objects_27.rds"))
dir.create(dirname(file.path(base, "output", "step26_figs_v2", "x")), showWarnings = FALSE)
saveRDS(list(BP = eBP, CC = eCC, MF = eMF, KEGG = eK), file.path(base, "output", "step26_figs_v2", "enrich_objects_27.rds"))
cat("STEP26 DONE\n")


############################################################
# FILE: step26_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step26_fig2_fig3_v2.R


############################################################
# FILE: step26b_cnet_v2.R
############################################################
# Step 26b: cnetplots only (fixed for current enrichplot API)
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(enrichplot); library(clusterProfiler); library(ggplot2) })
base <- "PROJECT_ROOT"
figdir <- file.path(base, "output", "paper", "figures_v2")
objs <- readRDS(file.path(base, "output", "step26_figs_v2", "enrich_objects_27.rds"))
grd27 <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
grd27 <- grd27[grd27$gene != "MIR644A", ]
library(org.Hs.eg.db); library(AnnotationDbi)
eg <- bitr(grd27$gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
eg <- eg[!duplicated(eg$SYMBOL), ]
fc_id <- grd27$logFC[match(eg$SYMBOL, grd27$gene)]; names(fc_id) <- eg$ENTREZID
for (nm in names(objs)) {
  e <- objs[[nm]]
  if (is.null(e) || nrow(as.data.frame(e)) == 0) { cat(nm, "empty, skip\n"); next }
  svg(file.path(figdir, paste0("Fig3_cnetplot_", nm, ".svg")), width = 9, height = 8)
  tryCatch({
    print(cnetplot(e, showCategory = 5, foldChange = fc_id, node_label = "all") +
          ggtitle(paste0("Top 5 ", nm, " terms (27 protein-coding GRDEGs)")))
    cat("saved cnetplot", nm, "\n")
  }, error = function(err) cat("cnetplot", nm, "failed:", conditionMessage(err), "\n"))
  dev.off()
}
cat("STEP26B DONE\n")


############################################################
# FILE: step27_fig4_fig5_v2.R
############################################################
# Step 27: Fig4A-D (27-gene selection) + Fig5A-F (lambda.min panel evaluation) -> figures_v2
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(glmnet); library(ggplot2); library(rms); library(pROC)
})
base <- "PROJECT_ROOT"
figdir <- file.path(base, "output", "paper", "figures_v2")
dir.create(figdir, recursive = TRUE, showWarnings = FALSE)
s13 <- file.path(base, "output", "step13_formal_rerun")

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"), check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
y <- factor(pheno$group, levels = c("Control", "DN"))

# ---- Fig4A forest (27 univariate ORs) ----
uni <- read.csv(file.path(s13, "S2_univariate_logistic_27GRDEG.csv"), stringsAsFactors = FALSE)
uni <- uni[order(uni$OR), ]
uni$gene <- factor(uni$gene, levels = uni$gene)
pA <- ggplot(uni, aes(x = OR, y = gene)) +
  geom_point(size = 2.2, color = "#2166AC") +
  geom_errorbarh(aes(xmin = pmax(CI_low, 1e-3), xmax = pmin(CI_high, 1e3)), height = 0.22, color = "#2166AC") +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey50") +
  scale_x_log10(limits = c(1e-3, 1e3)) +
  labs(title = "Univariate logistic regression (27 protein-coding GRDEGs, all p<0.05)",
       x = "Odds Ratio (log10)", y = NULL) + theme_bw(base_size = 10)
ggsave(file.path(figdir, "Fig4A_forest_univariate.svg"), pA, width = 8, height = 7, device = grDevices::svg)

# ---- Fig4B RF importance (27) ----
imp <- read.csv(file.path(s13, "RF_importance_27GRDEG.csv"), stringsAsFactors = FALSE)
imp$gene <- factor(imp$gene, levels = rev(imp$gene))
top10 <- head(imp[order(-imp$MeanDecreaseGini), ], 10)$gene
pB <- ggplot(imp, aes(x = MeanDecreaseGini, y = gene, fill = gene %in% top10)) +
  geom_col(width = 0.68) +
  scale_fill_manual(values = c(`TRUE` = "#B2182B", `FALSE` = "grey75"), guide = "none") +
  labs(title = "Random forest importance (mean decrease in Gini; top 10 highlighted)",
       x = "Mean decrease in Gini", y = NULL) + theme_bw(base_size = 10)
ggsave(file.path(figdir, "Fig4B_rf_importance.svg"), pB, width = 8, height = 7, device = grDevices::svg)

# ---- Fig4C LASSO CV curve + Fig4D path (from saved cv object) ----
cv <- readRDS(file.path(s13, "cv_glmnet.rds"))
svg(file.path(figdir, "Fig4C_lasso_cv.svg"), width = 7.5, height = 6)
plot(cv, main = "LASSO 10-fold cross-validation (27 GRDEGs)")
abline(v = log(c(cv$lambda.min, cv$lambda.1se)), lty = 2, col = c("#B2182B", "#2166AC"))
legend("topright", legend = c(sprintf("lambda.min=%.4f", cv$lambda.min), sprintf("lambda.1se=%.4f", cv$lambda.1se)),
       col = c("#B2182B", "#2166AC"), lty = 2, cex = 0.8)
dev.off()
svg(file.path(figdir, "Fig4D_lasso_trajectory.svg"), width = 7.5, height = 6)
plot(cv$glmnet.fit, xvar = "lambda", label = TRUE, main = "LASSO coefficient trajectory")
abline(v = log(cv$lambda.min), lty = 2)
dev.off()

# ---- Fig5A-F: lambda.min panel (FBN1, FN1, LPL, SLIT3, CCL21, SELL, EGF) ----
cf5 <- read.csv(file.path(s13, "LASSO_genes_lambda.min.csv"), stringsAsFactors = FALSE)
mg <- cf5$gene; cf <- cf5$coefficient; names(cf) <- mg
cat("lambda.min panel:", paste(mg, collapse = ","), "\n")
df <- as.data.frame(t(expr[mg, ])); df$y <- y
rs <- as.numeric(as.matrix(df[, mg]) %*% cf[mg])
df$RiskScore <- rs

# A nomogram
dd <- datadist(df); options(datadist = "dd")
f <- lrm(y ~ RiskScore, data = df, x = TRUE, y = TRUE)
svg(file.path(figdir, "Fig5A_nomogram_training.svg"), width = 11, height = 6)
plot(nomogram(f, fun = plogis, funlabel = "P(DN)", fun.at = c(0.01, 0.5, 0.99), nint = 5),
     main = "Nomogram: lambda.min 7-gene panel (GSE96804)", cex.axis = 0.75)
dev.off()
# B ROC
roc_obj <- roc(response = y, predictor = rs, quiet = TRUE)
ci <- ci.auc(roc_obj)
svg(file.path(figdir, "Fig5B_roc_training.svg"), width = 6.5, height = 6)
plot(roc_obj, col = "#B2182B", lwd = 2.2,
     main = sprintf("ROC: lambda.min 7-gene panel (AUC=%.4f, 95%% CI %.4f-%.4f)", auc(roc_obj), ci[1], ci[3]),
     print.auc = FALSE)
dev.off()
# C calibration (bootstrap 300)
set.seed(2024)
cal <- suppressWarnings(calibrate(f, method = "boot", B = 300))
svg(file.path(figdir, "Fig5C_calibration_training.svg"), width = 7, height = 7)
plot(cal, main = "Bootstrap calibration (B = 300), lambda.min panel")
dev.off()
# D DCA
pr <- plogis(rs)
prev <- mean(y == "DN")
thr <- seq(0.01, 0.99, by = 0.01)
nb_model <- sapply(thr, function(t) {
  tp <- mean(pr >= t & y == "DN"); fp <- mean(pr >= t & y == "Control")
  tp - fp * (t / (1 - t)) })
nb_all <- sapply(thr, function(t) prev - (1 - prev) * (t / (1 - t)))
dca <- data.frame(threshold = rep(thr, 3),
                  NB = c(nb_model, nb_all, rep(0, length(thr))),
                  Strategy = rep(c("Model (lambda.min 7-gene)", "Treat all", "Treat none"), each = length(thr)))
pD <- ggplot(dca, aes(x = threshold, y = NB, color = Strategy)) +
  geom_line(linewidth = 0.9) + coord_cartesian(ylim = c(-0.05, max(nb_model) + 0.05)) +
  labs(title = "Decision curve analysis (training set)", x = "Threshold probability", y = "Net benefit") +
  theme_bw(base_size = 10) + theme(legend.position = "bottom")
ggsave(file.path(figdir, "Fig5D_dca_training.svg"), pD, width = 7, height = 5, device = grDevices::svg)
# E friends
fr <- read.csv(file.path(base, "output", "step16_downstream_v2", "friends_mean_similarity_final7.csv"), stringsAsFactors = FALSE)
fr <- fr[order(-fr$mean_similarity), ]
fr$gene <- factor(fr$gene, levels = rev(fr$gene))
pE <- ggplot(fr, aes(x = mean_similarity, y = gene)) +
  geom_col(fill = "#4393C3", width = 0.65) +
  labs(title = "GOSemSim friends analysis (mean GO-BP semantic similarity)",
       subtitle = "FN1 is the most functionally coherent member", x = "Mean semantic similarity", y = NULL) +
  theme_bw(base_size = 10)
ggsave(file.path(figdir, "Fig5E_gosemsim_friends.svg"), pE, width = 7, height = 6, device = grDevices::svg)
# F risk score box
dfb <- data.frame(RiskScore = rs, Group = y)
pF <- ggplot(dfb, aes(x = Group, y = RiskScore, fill = Group)) +
  geom_boxplot(width = 0.45, outlier.shape = NA, alpha = 0.85) +
  geom_jitter(width = 0.12, size = 1.1, alpha = 0.6) +
  scale_fill_manual(values = c(Control = "#4393C3", DN = "#D6604D"), guide = "none") +
  labs(title = "RiskScore by group (lambda.min 7-gene panel)",
       subtitle = sprintf("Mann-Whitney p = 3.5e-10; mean DN vs Control = %.2f vs %.2f",
                          mean(rs[y == "DN"]), mean(rs[y == "Control"])),
       x = NULL, y = "RiskScore") + theme_bw(base_size = 10)
ggsave(file.path(figdir, "Fig5F_riskscore_box_training.svg"), pF, width = 6.5, height = 5.5, device = grDevices::svg)
write.csv(data.frame(gsm = pheno$gsm, group = pheno$group, RiskScore = rs),
          file.path(base, "output", "step27_figs_v2_riskscore_training_lambda_min.csv"), row.names = FALSE)
cat("STEP27 DONE\n")


############################################################
# FILE: step28_fig678_v2.R
############################################################
# Step 28: Fig6/7 (external validation, lambda.min frozen coefficients) + Fig8A-D -> figures_v2
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(rms); library(pROC); library(ggplot2); library(pheatmap) })
base <- "PROJECT_ROOT"
figdir <- file.path(base, "output", "paper", "figures_v2")
s13 <- file.path(base, "output", "step13_formal_rerun")

cf <- read.csv(file.path(s13, "LASSO_genes_lambda.min.csv"), stringsAsFactors = FALSE)
mg <- cf$gene; coefv <- cf$coefficient; names(coefv) <- mg
expr_tr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"), check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr_tr) <- expr_tr$gene; expr_tr$gene <- NULL; expr_tr <- as.matrix(expr_tr)
ph_tr <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
ytr <- factor(ph_tr$group, levels = c("Control", "DN"))

mk_eval <- function(gse, tag, exprm, phv) {
  yv <- factor(phv$group, levels = c("Control", if (gse == "GSE30528") "DKD" else "DKD"))
  ev <- as.data.frame(t(exprm[mg, colnames(exprm)])); ev$y <- yv
  rs <- as.numeric(as.matrix(ev[, mg]) %*% coefv[mg])
  ev$RiskScore <- rs
  dd <- datadist(ev); options(datadist = dd)   # pass object directly (local scope)
  f <- suppressWarnings(lrm(y ~ RiskScore, data = ev, x = TRUE, y = TRUE))
  # A nomogram
  svg(file.path(figdir, paste0("Fig", tag, "A_nomogram_", gse, ".svg")), width = 11, height = 6)
  plot(nomogram(f, fun = plogis, funlabel = "P(DKD)", fun.at = c(0.01, 0.5, 0.99), nint = 5),
       main = paste0("Nomogram: lambda.min 7-gene panel (", gse, ")"), cex.axis = 0.75)
  dev.off()
  # B ROC
  r <- roc(response = yv, predictor = rs, quiet = TRUE); ci <- ci.auc(r)
  mt <- suppressWarnings(wilcox.test(rs ~ yv, exact = FALSE)$p.value)
  svg(file.path(figdir, paste0("Fig", tag, "B_roc_", gse, ".svg")), width = 6.5, height = 6)
  plot(r, col = "#B2182B", lwd = 2.2,
       main = sprintf("%s: AUC=%.3f (95%% CI %.3f-%.3f), MWU p=%.3f", gse, auc(r), ci[1], ci[3], mt))
  dev.off()
  # C calibration
  set.seed(2024)
  cal <- suppressWarnings(calibrate(f, method = "boot", B = 200))
  svg(file.path(figdir, paste0("Fig", tag, "C_calibration_", gse, ".svg")), width = 7, height = 7)
  plot(cal, main = paste0("Bootstrap calibration (B=200), ", gse))
  dev.off()
  # D DCA
  pr <- plogis(rs); prev <- mean(yv == "DKD")
  thr <- seq(0.01, 0.99, by = 0.01)
  nb_m <- sapply(thr, function(t) mean(pr >= t & yv == "DKD") - mean(pr >= t & yv == "Control") * (t / (1 - t)))
  nb_a <- sapply(thr, function(t) prev - (1 - prev) * (t / (1 - t)))
  dca <- data.frame(threshold = rep(thr, 3), NB = c(nb_m, nb_a, rep(0, length(thr))),
                    Strategy = rep(c("Model (lambda.min 7-gene)", "Treat all", "Treat none"), each = length(thr)))
  pD <- ggplot(dca, aes(x = threshold, y = NB, color = Strategy)) +
    geom_line(linewidth = 0.9) +
    coord_cartesian(ylim = c(-0.05, max(nb_m, nb_a) + 0.05)) +
    labs(title = paste0("Decision curve analysis (", gse, ")"), x = "Threshold probability", y = "Net benefit") +
    theme_bw(base_size = 10) + theme(legend.position = "bottom")
  ggsave(file.path(figdir, paste0("Fig", tag, "D_dca_", gse, ".svg")), pD, width = 7, height = 5, device = grDevices::svg)
  cat(gse, "AUC", round(auc(r), 3), "p", signif(mt, 3), "\n")
  list(rs = rs, y = yv)
}

# ---- Fig6: GSE30528 / Fig7: GSE30529 ----
for (gse in c("GSE30528", "GSE30529")) {
  ev <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")), check.names = FALSE, stringsAsFactors = FALSE)
  rownames(ev) <- ev$gene; ev$gene <- NULL; ev <- as.matrix(ev)
  phv <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  mk_eval(gse, if (gse == "GSE30528") 6 else 7, ev, phv)
}

# ---- Fig8A single-gene ROC (training, 7 genes) ----
cols <- rainbow(length(mg)); ltyv <- rep(1:2, length.out = length(mg))
svg(file.path(figdir, "Fig8A_single_gene_roc_training.svg"), width = 6.5, height = 6)
plot(x = NULL, y = NULL, xlim = c(1, 0), ylim = c(0, 1),
     xlab = "1 - Specificity", ylab = "Sensitivity",
     main = "Single-gene ROC curves (training set GSE96804, 7-gene panel)", cex.main = 0.95)
leg_lab <- character(0); leg_col <- character(0); leg_lty <- numeric(0)
for (i in seq_along(mg)) {
  r <- roc(response = ytr, predictor = expr_tr[mg[i], ], quiet = TRUE)
  lines(r, col = cols[i], lwd = 1.8, lty = ltyv[i])
  leg_lab <- c(leg_lab, sprintf("%s (AUC=%.3f)", mg[i], auc(r)))
  leg_col <- c(leg_col, cols[i]); leg_lty <- c(leg_lty, ltyv[i])
}
legend("bottomright", legend = leg_lab, col = leg_col, lwd = 1.8, lty = leg_lty, cex = 0.72, bty = "n")
dev.off()

# ---- Fig8B correlation matrix (lambda.min 7 genes) ----
cm <- suppressWarnings(cor(t(expr_tr[mg, ]), method = "spearman"))
svg(file.path(figdir, "Fig8B_correlation_matrix.svg"), width = 7, height = 6)
pheatmap(cm, color = colorRampPalette(c("#2166AC", "white", "#B2182B"))(100),
         display_numbers = TRUE, number_size = 7, fontsize = 9, border_color = "grey80",
         main = "Spearman correlations among panel genes (FBN1-FN1 r=0.899)")
dev.off()

# ---- Fig8C/D single-gene expression boxes in validation cohorts ----
for (gse in c("GSE30528", "GSE30529")) {
  ev <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")), check.names = FALSE, stringsAsFactors = FALSE)
  rownames(ev) <- ev$gene; ev$gene <- NULL
  phv <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  yv <- factor(phv$group, levels = c("Control", "DKD"))
  dfb <- data.frame(gene = rep(mg, each = length(yv)),
                    expr = as.numeric(t(ev[mg, ])),
                    Group = rep(yv, times = length(mg)))
  dfb$gene <- factor(dfb$gene, levels = mg)
  p8 <- ggplot(dfb, aes(x = gene, y = expr, fill = Group)) +
    geom_boxplot(width = 0.55, outlier.shape = NA, alpha = 0.85) +
    geom_jitter(width = 0.1, size = 0.7, alpha = 0.55) +
    scale_fill_manual(values = c(Control = "#4393C3", DKD = "#D6604D")) +
    labs(title = paste0("Panel-gene expression by group (", gse, ")"),
         subtitle = "Directionally consistent: LPL/FBN1/EGF/FN1; heterogeneous: SELL/SLIT3/CCL21",
         x = NULL, y = "Normalized expression (centered)") +
    theme_bw(base_size = 10) + theme(axis.text.x = element_text(angle = 30, hjust = 1), legend.position = "bottom")
  ggsave(file.path(figdir, paste0("Fig8", if (gse == "GSE30528") "C" else "D", "_modelgene_box_", gse, ".svg")),
         p8, width = 8, height = 5.5, device = grDevices::svg)
}
cat("STEP28 DONE\n")


############################################################
# FILE: step29_fig9_svr.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 29: Fig9 SVR-primary immune figures (A stacked, B cell-corr heatmap,
C gene-cell bubble) + ssGSEA comparator renamed for supplementary. -> figures_v2"""
import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.stats import spearmanr

BASE = "PROJECT_ROOT"
FIG = os.path.join(BASE, "output", "paper", "figures_v2")
os.makedirs(FIG, exist_ok=True)

prop = pd.read_csv(os.path.join(BASE, "output/step17_cibersort_svr/cibersort_proportions_official_LM22.csv"),
                   index_col=0, sep=None, engine="python")
pheno = pd.read_csv(os.path.join(BASE, "data/GSE96804_pheno.csv"))
groups = pheno.set_index("gsm")["group"]
prop = prop.loc[[g for g in pheno["gsm"] if g in prop.index]]
grp = groups.loc[prop.index]
corr = prop.corr(method="spearman")
panel = ["FBN1", "FN1", "LPL", "SLIT3", "CCL21", "SELL", "EGF"]
expr = pd.read_csv(os.path.join(BASE, "data/GSE96804_norm_expr.csv"), index_col=0)
expr.columns = [c.strip() for c in expr.columns]
common = [c for c in prop.index if c in expr.columns]
gc = pd.read_csv(os.path.join(BASE, "output/step17_cibersort_svr/gene_celltype_SVR_corr_full.csv"))
gc = gc[gc.gene.isin(panel)]

# ---- Fig9A stacked bars of SVR fractions ----
order = prop.sum().sort_values(ascending=False).index
top12 = list(order[:12]); rest = [c for c in prop.columns if c not in top12]
plotdf = prop[top12].copy()
plotdf["Other (10 types)"] = prop[rest].sum(axis=1)
cols = plt.cm.tab20(np.linspace(0, 1, plotdf.shape[1]))
fig, ax = plt.subplots(figsize=(11, 5.5))
bottom = np.zeros(len(plotdf))
sample_color = {"DN": "#C0392B", "Control": "#2166AC"}
for i, c in enumerate(plotdf.columns):
    ax.bar(range(len(plotdf)), plotdf[c], bottom=bottom, color=cols[i], width=0.82, label=c)
    bottom += plotdf[c].values
ax.set_xticks(range(len(plotdf)))
ax.set_xticklabels([f"{s[:10]}\n{groups[s]}" for s in plotdf.index], fontsize=6, rotation=90)
ax.set_ylabel("Deconvolved fraction (nu-SVR, official LM22)")
ax.set_title("Immune-cell proportions per sample (GSE96804) - CIBERSORT nu-SVR, zero-abundance types retained as 0")
ax.legend(fontsize=6.5, ncol=3, loc="upper center", bbox_to_anchor=(0.5, -0.28))
plt.tight_layout()
plt.savefig(os.path.join(FIG, "Fig9A_immune_stacked_SVR.svg"), bbox_inches="tight")
plt.close()

# ---- Fig9B cell-type correlation heatmap (SVR) ----
order2 = corr.mean().sort_values(ascending=False).index
cm = corr.loc[order2, order2]
fig, ax = plt.subplots(figsize=(9.5, 8))
im = ax.imshow(cm.values, cmap="RdBu_r", vmin=-0.7, vmax=0.7)
ax.set_xticks(range(len(cm))); ax.set_yticks(range(len(cm)))
ax.set_xticklabels(cm.columns, rotation=90, fontsize=7)
ax.set_yticklabels(cm.index, fontsize=7)
for i in range(len(cm)):
    for j in range(len(cm)):
        ax.text(j, i, f"{cm.values[i, j]:.2f}", ha="center", va="center", fontsize=5.2,
                color="white" if abs(cm.values[i, j]) > 0.45 else "black")
ax.set_title("Spearman correlations among 22 cell-type fractions (nu-SVR)\n"
             "96 positive / 135 negative of 231 pairs - rank-method all-positive artifact abolished", fontsize=9.5)
plt.colorbar(im, shrink=0.75)
plt.tight_layout()
plt.savefig(os.path.join(FIG, "Fig9B_celltype_corr_heatmap_SVR.svg"), bbox_inches="tight")
plt.close()

# ---- Fig9C gene-cell bubble (SVR, BH-corrected) ----
gc["absr"] = gc.spearman_r.abs()
gc = gc.sort_values("absr", ascending=False)
cell_order = (gc.groupby("cell_type").spearman_r.apply(lambda s: s.abs().max())
                .sort_values(ascending=False).head(14).index.tolist())
gc2 = gc[gc.cell_type.isin(cell_order)].copy()
piv = gc2.pivot(index="gene", columns="cell_type", values="spearman_r").reindex(panel)[cell_order]
size = gc2.pivot(index="gene", columns="cell_type", values="absr").reindex(panel)[cell_order] * 260
fig, ax = plt.subplots(figsize=(10.5, 4.6))
for gi, g in enumerate(piv.index):
    for ci, c in enumerate(piv.columns):
        v = piv.loc[g, c]
        if np.isnan(v):
            continue
        ax.scatter(ci, gi, s=size.loc[g, c], c="#B2182B" if v > 0 else "#2166AC",
                   alpha=0.85, edgecolors="none")
        ax.text(ci, gi, f"{v:.2f}", ha="center", va="center", fontsize=5.6,
                color="white" if abs(v) > 0.55 else "black")
ax.set_yticks(range(len(piv.index))); ax.set_yticklabels(piv.index, fontsize=8)
ax.set_xticks(range(len(piv.columns)))
ax.set_xticklabels([c.replace(".", " ") for c in piv.columns], rotation=35, ha="right", fontsize=7.5)
ax.axhline(len(panel) - 1 - list(piv.index).index("SELL") - 0.0, color="none")
ax.set_title("Panel gene x cell-type Spearman r (nu-SVR fractions)\n"
             "SELL-neutrophil r=0.308 (0.327 decontaminated); 113/154 pairs q<0.05 (BH)", fontsize=9.5)
plt.tight_layout()
plt.savefig(os.path.join(FIG, "Fig9C_gene_celltype_bubble_SVR.svg"), bbox_inches="tight")
plt.close()
print("Fig9 SVR versions saved")

# ---- Supplementary S9B-ssGSEA heatmap (rank-method comparator) ----
try:
    nes = pd.read_csv(os.path.join(BASE, "output/step9/immune_nes_matrix.csv"), index_col=0)
    nes.columns = [c.strip() for c in nes.columns]
    cnes = nes.corr(method="spearman")
    o3 = cnes.mean().sort_values(ascending=False).index
    cnes = cnes.loc[o3, o3]
    fig, ax = plt.subplots(figsize=(9.5, 8))
    im = ax.imshow(cnes.values, cmap="RdBu_r", vmin=0, vmax=1)
    ax.set_xticks(range(len(cnes))); ax.set_yticks(range(len(cnes)))
    ax.set_xticklabels(cnes.columns, rotation=90, fontsize=7)
    ax.set_yticklabels(cnes.index, fontsize=7)
    ax.set_title("SUPPLEMENTARY comparator: ssGSEA NES correlations (all positive)\n"
                 "uniformly positive structure = rank-method artifact (superseded)", fontsize=9.5)
    plt.colorbar(im, shrink=0.75)
    plt.tight_layout()
    plt.savefig(os.path.join(FIG, "SuppFigS9B_celltype_corr_heatmap_ssGSEA_comparator.svg"), bbox_inches="tight")
    plt.close()
    print("ssGSEA comparator heatmap saved")
except Exception as e:
    print("ssGSEA comparator failed:", e)
print("STEP29 DONE")


############################################################
# FILE: step2a_uniprot_grg.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 2a: Fetch human reviewed glycosylation-related genes from UniProt REST.
Keyword KW-0325 = Glycosylation (mechanism/enzymes). Save gene symbols with accessions.
"""
import json, os, sys, time
import urllib.request, urllib.parse, ssl

CTX = ssl.create_default_context()
CTX.check_hostname = False
CTX.verify_mode = ssl.CERT_NONE

BASE = "https://rest.uniprot.org/uniprotkb/search"
QUERY = "(keyword:KW-0325) AND organism_id:9606 AND reviewed:true"
FIELDS = "accession,id,gene_primary,protein_name"

def fetch(url, tries=4):
    last = None
    for i in range(tries):
        try:
            req = urllib.request.Request(url, headers={"User-Agent": "dn-glyco-repro/1.0"})
            with urllib.request.urlopen(req, timeout=90, context=CTX) as r:
                return r.read().decode("utf-8")
        except Exception as e:
            last = e
            time.sleep(2 * (i + 1))
    raise RuntimeError(f"fetch failed: {url} :: {last}")

def main():
    out_dir = "PROJECT_ROOT/data"
    os.makedirs(out_dir, exist_ok=True)
    cursor = None
    rows = []
    page = 0
    while True:
        params = {
            "query": QUERY,
            "fields": FIELDS,
            "format": "tsv",
            "size": "500",
        }
        if cursor:
            params["cursor"] = cursor
        url = BASE + "?" + urllib.parse.urlencode(params)
        text = fetch(url)
        lines = [l for l in text.strip().split("\n") if l.strip()]
        if page == 0 and lines:
            header = lines[0]
            rows.append(("__header__", header))
        for l in lines[1:]:
            parts = l.split("\t")
            acc = parts[0] if len(parts) > 0 else ""
            entry = parts[1] if len(parts) > 1 else ""
            gene = parts[2] if len(parts) > 2 else ""
            pname = parts[3] if len(parts) > 3 else ""
            rows.append((gene, acc + "\t" + entry + "\t" + pname))
        page += 1
        # parse Link header for next cursor via response headers is easier with urllib:
        # do a second request trick: UniProt returns Link rel=next; urllib exposes headers only during read.
        # Simpler: request with cursor loop until fewer than size rows returned.
        n_rows = len(lines) - 1
        print(f"page {page}: {n_rows} rows", flush=True)
        if n_rows < 500:
            break
        # get next cursor: need Link header; redo request capturing headers
        req = urllib.request.Request(url, headers={"User-Agent": "dn-glyco-repro/1.0"})
        with urllib.request.urlopen(req, timeout=90, context=CTX) as r:
            link = r.headers.get("Link", "")
        nxt = None
        for part in link.split(","):
            if 'rel="next"' in part:
                nxt = part.split("<")[1].split(">")[0]
        if not nxt:
            break
        cursor = nxt.split("cursor=")[-1]
        time.sleep(1)
    out_txt = os.path.join(out_dir, "uniprot_kw0325_human.tsv")
    with open(out_txt, "w") as f:
        for gene, rest in rows:
            f.write(gene + "\t" + rest + "\n")
    genes = set()
    for gene, _ in rows:
        if gene and gene != "__header__":
            for g in gene.replace(" ", "").split(";"):
                if g:
                    genes.add(g)
    print("unique gene symbols:", len(genes))
    with open(os.path.join(out_dir, "uniprot_kw0325_genes.txt"), "w") as f:
        f.write("\n".join(sorted(genes)) + "\n")
    print("saved:", out_txt)

if __name__ == "__main__":
    main()


############################################################
# FILE: step2b_go_grg.R
############################################################
# Step 2b: Build GO-based glycosylation-related gene set (terms + all offspring)
# Run with: /public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step2b_go_grg.R
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(GO.db); library(org.Hs.eg.db); library(AnnotationDbi)
})
out_dir <- "PROJECT_ROOT/data"

# Core glycosylation-related GO terms (mirroring GeneCards "Glycosylation" semantics)
bp_ids <- c("GO:0006486",  # protein N-linked glycosylation
            "GO:0006493",  # protein O-linked glycosylation
            "GO:0009101",  # glycoprotein biosynthetic process
            "GO:0006022")  # aminoglycan (glycosaminoglycan) metabolic process
mf_ids <- c("GO:0016757",  # glycosyltransferase activity
            "GO:0016798",  # glycosidase activity
            "GO:0005539")  # glycosaminoglycan binding

bp_off <- unique(unlist(as.list(GOBPOFFSPRING[bp_ids])))
mf_off <- unique(unlist(as.list(GOMFOFFSPRING[mf_ids])))
bp_all <- unique(c(bp_ids, bp_off)); bp_all <- bp_all[!is.na(bp_all)]
mf_all <- unique(c(mf_ids, mf_off)); mf_all <- mf_all[!is.na(mf_all)]
cat("BP terms incl offspring:", length(bp_all), " MF terms incl offspring:", length(mf_all), "\n")

g2a <- toTable(org.Hs.egGO2ALLEGS)
sel <- g2a[(g2a$go_id %in% bp_all & g2a$Ontology == "BP") |
           (g2a$go_id %in% mf_all & g2a$Ontology == "MF"), ]
cat("annotations rows:", nrow(sel), "\n")
sym_map <- toTable(org.Hs.egSYMBOL)
m <- merge(sel[, c("gene_id", "go_id", "Ontology")], sym_map, by.x = "gene_id", by.y = "gene_id", all.x = TRUE)
m <- m[!is.na(m$symbol) & m$symbol != "", ]
m <- m[!duplicated(paste(m$symbol, m$go_id)), ]
res <- m[, c("symbol", "go_id", "Ontology")]
res <- res[order(res$symbol), ]
write.table(res, file.path(out_dir, "go_grg_annotations.tsv"),
            sep = "\t", quote = FALSE, row.names = FALSE)
genes <- sort(unique(res$symbol))
cat("unique GO-GRG gene symbols:", length(genes), "\n")
writeLines(genes, file.path(out_dir, "go_grg_genes.txt"))
cat("saved go_grg_annotations.tsv / go_grg_genes.txt\n")


############################################################
# FILE: step30_additional_files.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Assemble Additional files 1-16: compute missing pieces (AF3, AF5, AF11),
merge AF8, copy all to output/additional_files/ and zip."""
import os
import numpy as np
import pandas as pd
from scipy.stats import spearmanr

BASE = "PROJECT_ROOT"
OUT = os.path.join(BASE, "output", "additional_files")
os.makedirs(OUT, exist_ok=True)

def bh(p):
    p = np.asarray(p, dtype=float); n = len(p)
    o = np.argsort(p); ranked = p[o] * n / np.arange(1, n + 1)
    q = np.minimum.accumulate(ranked[::-1])[::-1]
    out = np.empty(n); out[o] = np.clip(q, 0, 1)
    return out

panel = ["FBN1", "FN1", "LPL", "SLIT3", "CCL21", "SELL", "EGF"]

# ---- AF11: SVR gene-cell 154 pairs with p and BH q ----
gc = pd.read_csv(os.path.join(BASE, "output/step17_cibersort_svr/gene_celltype_SVR_corr_full.csv"))
expr = pd.read_csv(os.path.join(BASE, "data/GSE96804_norm_expr.csv"), index_col=0)
expr.columns = [c.strip() for c in expr.columns]
prop = pd.read_csv(os.path.join(BASE, "output/step17_cibersort_svr/cibersort_proportions_official_LM22.csv"),
                   index_col=0, sep=None, engine="python")
common = [c for c in expr.columns if c in prop.index]
rows = []
for g in panel:
    for ct in prop.columns:
        r, p = spearmanr(expr.loc[g, common].values.astype(float), prop.loc[common, ct].values.astype(float))
        rows.append(dict(gene=g, cell_type=ct, spearman_r=r, p_value=p))
af11 = pd.DataFrame(rows)
af11["q_BH"] = bh(af11.p_value.values)
af11.to_csv(os.path.join(OUT, "Additional_file_11_gene_celltype_nuSVR_corr_154pairs_BH.csv"), index=False)
print("AF11 done:", len(af11), "pairs; significant q<0.05:", (af11.q_BH < 0.05).sum())

# ---- AF3: ssGSEA + nu-SVR gene-cell matrices + ssGSEA cell-cell (Fig9b comparator) ----
tmp3 = os.path.join(OUT, "Additional_file_3_tmp"); os.makedirs(tmp3, exist_ok=True)
af11.to_csv(os.path.join(tmp3, "AF3a_model_gene_celltype_nuSVR_corr_154pairs_BH.csv"), index=False)
s3g = pd.read_csv(os.path.join(BASE, "output/../supplementary/S3_gene_immune_corr_BH_newmodel.csv"))
s3g.to_csv(os.path.join(tmp3, "AF3b_model_gene_celltype_ssGSEA_corr_154pairs_BH.csv"), index=False)
print("AF3b ssGSEA pairs:", len(s3g))
nes = pd.read_csv(os.path.join(BASE, "output/step9/immune_nes_matrix.csv"), index_col=0)
nes.columns = [c.strip() for c in nes.columns]
cnes = nes.corr(method="spearman")
cnes.to_csv(os.path.join(tmp3, "AF3c_ssGSEA_celltype_corr_22x22_Fig9b_comparator.csv"))
svrcell = prop.corr(method="spearman")
svrcell.to_csv(os.path.join(tmp3, "AF3d_nuSVR_celltype_corr_22x22_Fig9b.csv"))
zf = os.path.join(OUT, "Additional_file_3_gene_immune_corr_matrices.zip")
if os.path.exists(zf): os.remove(zf)
os.system(f'cd "{tmp3}" && zip -q "{zf}" *.csv && cd .. && rm -rf "{tmp3}"')
print("AF3 zip done:", zf)

# ---- AF5: DCA net-benefit tables recomputed for lambda.min panel ----
cf = pd.read_csv(os.path.join(BASE, "output/step13_formal_rerun/LASSO_genes_lambda.min.csv"))
mg = cf.gene.tolist(); coefv = cf.coefficient.values; coefv = dict(zip(mg, coefv))
def rs_of(mat):
    return sum(coefv[g] * mat.loc[g].values for g in mg)
def dca_table(rs, y):
    pr = 1 / (1 + np.exp(-rs)); prev = np.mean(y == "DN")
    thr = np.arange(0.01, 1.0, 0.01)
    nb_m = [np.mean(pr >= t) * 0 for t in thr]  # placeholder replaced below
    nb_m = [np.mean((pr >= t) & (y == "DN")) - np.mean((pr >= t) & (y == "Control")) * (t / (1 - t)) for t in thr]
    nb_a = [prev - (1 - prev) * (t / (1 - t)) for t in thr]
    return pd.DataFrame(dict(threshold=thr, NB_model_lambda_min_7gene=nb_m, NB_treat_all=nb_a, NB_treat_none=0.0))

ph = pd.read_csv(os.path.join(BASE, "data/GSE96804_pheno.csv"))
rs_tr = pd.read_csv(os.path.join(BASE, "output/step27_figs_v2_riskscore_training_lambda_min.csv"))
d1 = dca_table(rs_tr.RiskScore.values, rs_tr.group.values)
d1.to_csv(os.path.join(OUT, "Additional_file_5a_DCA_netbenefit_training_lambda_min.csv"), index=False)
coefs = np.array([coefv[g] for g in mg])
for gse, tag in [("GSE30528", "5b"), ("GSE30529", "5c")]:
    ev = pd.read_csv(os.path.join(BASE, f"data/{gse}_norm_expr.csv"), index_col=0)
    phv = pd.read_csv(os.path.join(BASE, f"data/{gse}_pheno.csv"))
    yv = (phv.group == "DKD").astype(int).values
    ctl = (phv.group == "Control").values
    rs = coefs @ ev.loc[mg, ev.columns].values.astype(float)
    pr = 1 / (1 + np.exp(-rs)); prev = yv.mean()
    thr = np.arange(0.01, 1.0, 0.01)
    nb_m = [np.mean((pr >= t) & (yv == 1)) - np.mean((pr >= t) & (ctl)) * (t / (1 - t)) for t in thr]
    nb_a = [prev - (1 - prev) * (t / (1 - t)) for t in thr]
    pd.DataFrame(dict(threshold=thr, NB_model_lambda_min_7gene=nb_m, NB_treat_all=nb_a, NB_treat_none=0.0)) \
      .to_csv(os.path.join(OUT, f"Additional_file_{tag}_DCA_netbenefit_{gse}_lambda_min.csv"), index=False)
print("AF5 done (3 DCA tables, lambda.min panel)")

# ---- AF8: 231 pairs + partial correlations merged ----
cc = pd.read_csv(os.path.join(BASE, "output/step18_global_factor/cellcell_231pairs_BH.csv"))
pc = pd.read_csv(os.path.join(BASE, "output/step18_global_factor/cellcell_partialcorr_given_PC1.csv"))
print("AF8 partial cols:", list(pc.columns))
key = lambda df: df.cellA.astype(str) + "|" + df.cellB.astype(str)
if {"cellA", "cellB"}.issubset(pc.columns):
    m = cc.merge(pc, on=["cellA", "cellB"], how="left", suffixes=("", "_pc"))
else:  # long-format partial table (cellA|cellB keys inside)
    print("partial table head:", pc.head(2).to_dict("records"))
    m = cc
m.to_csv(os.path.join(OUT, "Additional_file_8_cellcell_231pairs_BH_with_partialcorr_PC1.csv"), index=False)
print("AF8 merged:", m.shape)

# ---- copy the rest ----
copies = {
 1: ("output/step12_coding_filter/GRG_list_v2_protein_coding.csv", "Additional_file_1_protein_coding_GRG_list_812.csv"),
 2: ("output/step13_formal_rerun/S2_univariate_logistic_27GRDEG.csv", "Additional_file_2_univariate_logistic_27GRDEGs.csv"),
 4: ("output/step16_downstream_v2/single_gene_AUC_final7.csv", "Additional_file_4_single_gene_AUC_final_panel.csv"),
 6: ("output/step21_leakfree_cv/leakfree_cv_oof_auc.csv", "Additional_file_6_leakfree_repeated_CV_OOF_AUCs.csv"),
 7: ("output/step14_robustness/bootstrap_selection_frequency_1000.csv", "Additional_file_7_bootstrap_selection_frequencies_1000.csv"),
 9: ("output/step20_1se_panel/AUC_1se5gene_panel.csv", "Additional_file_9_external_AUC_lambda1se_5gene.csv"),
 10: ("output/step20_1se_panel/model_comparison_v4_main1se_vs_sensitivity.csv", "Additional_file_10_dual_panel_comparison_4layers.csv"),
 12: ("output/step22_svr_audit/celltype_zero_abundance_audit.csv", "Additional_file_12_celltype_zero_abundance_audit.csv"),
 13: ("output/step25_concordance_v2/three_method_concordance.csv", "Additional_file_13_three_method_concordance_SELL_neutrophil.csv"),
 14: ("output/step25_concordance_v2/panel_tradeoff_min_vs_1se.csv", "Additional_file_14_panel_tradeoff_min_vs_1se.csv"),
 15: ("output/step24_gene_crosscohort/gene_direction_consistency.csv", "Additional_file_15_crosscohort_gene_direction_consistency.csv"),
 16: ("output/step23_svr_crosscohort/svr_sell_neutrophil_three_cohorts.csv", "Additional_file_16_three_cohort_nuSVR_SELL_neutrophil.csv"),
}
for n, (src, dst) in copies.items():
    pd.read_csv(os.path.join(BASE, src)).to_csv(os.path.join(OUT, dst), index=False)
print("copied", len(copies), "static files")

os.system(f'cd "{OUT}" && rm -f Additional_files_1-16.zip && zip -q Additional_files_1-16.zip *.csv Additional_file_3*.zip && ls')
print("DONE")


############################################################
# FILE: step31_fig1_v2.R
############################################################
# Step 31: Regenerate Fig1 (normalization boxplots) without label overlap -> figures_v2
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(limma) })
avereps2 <- function(x) {
  out <- aggregate(x, by = list(sym = rownames(x)), FUN = mean)
  rn <- out$sym; out$sym <- NULL
  rownames(out) <- rn; as.matrix(out)
}
base <- "PROJECT_ROOT"
figdir <- file.path(base, "output", "paper", "figures_v2")
dir.create(figdir, recursive = TRUE, showWarnings = FALSE)

# ---- annotation (same as step3b) ----
ann571 <- read.delim(file.path(base, "data", "GPL571_table.txt"), skip = 1,
                     sep = "\t", quote = "\"", check.names = FALSE, stringsAsFactors = FALSE)
ann571 <- ann571[ann571[[1]] != "!platform_table_end", ]
g571 <- data.frame(probe = ann571$ID,
                   symbol = sub("\\s*$", "", sub("^\\s*", "", sapply(strsplit(as.character(ann571[["Gene Symbol"]]), "///"), function(x) trimws(x[1])))),
                   stringsAsFactors = FALSE)
ann17586 <- read.delim(file.path(base, "data", "GPL17586_table.txt"), skip = 1,
                       sep = "\t", quote = "\"", check.names = FALSE, stringsAsFactors = FALSE)
ann17586 <- ann17586[ann17586[[1]] != "!platform_table_end", ]
ga <- as.character(ann17586[["gene_assignment"]])
first_entry <- sub("\\s*///.*$", "", ga)
sym17586 <- sub("^[^/]*//\\s*", "", first_entry)
sym17586 <- trimws(sub("//.*$", "", paste0(sym17586, "//")))
g17586 <- data.frame(probe = ann17586$ID, symbol = trimws(sym17586), stringsAsFactors = FALSE)

plats <- list(GSE96804 = g17586, GSE30528 = g571, GSE30529 = g571)
boxdata <- list()
for (gse in names(plats)) {
  expr <- read.delim(file.path(base, "data", paste0(gse, "_expr.tsv")),
                     sep = "\t", header = TRUE, check.names = FALSE, stringsAsFactors = FALSE)
  rownames(expr) <- expr[[1]]; expr[[1]] <- NULL
  expr <- as.matrix(expr); storage.mode(expr) <- "numeric"
  pheno <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  stopifnot(identical(as.character(pheno$gsm), colnames(expr)))
  ann <- plats[[gse]]
  m <- merge(data.frame(probe = rownames(expr), stringsAsFactors = FALSE), ann, by = "probe", sort = FALSE)
  em <- expr[match(m$probe, rownames(expr)), , drop = FALSE]
  m$symbol <- ifelse(is.na(m$symbol) | m$symbol == "", paste0(m$probe, "_NA"), m$symbol)
  rownames(em) <- m$symbol
  em <- avereps2(em)
  if (max(em, na.rm = TRUE) > 100) em <- log2(em + 1)
  boxdata[[gse]] <- list(before = em, after = normalizeBetweenArrays(em), pheno = pheno)
  cat(gse, "ready:", dim(em), "\n")
}

draw_pair <- function(gse, bd) {
  ph <- bd$pheno
  ord <- order(ph$group)
  grp <- ph$group[ord]
  cols <- ifelse(grp %in% c("DN", "DKD"), "#E69F00", "#56B4E9")
  lbl_case <- if (any(grp == "DN")) "DN" else "DKD"
  n_case <- sum(grp %in% c("DN", "DKD")); n_ctl <- sum(grp == "Control")
  for (what in c("before", "after")) {
    em <- bd[[what]][, ord]
    boxplot(em, outline = FALSE, col = cols, pch = ".", names = NA,
            xaxt = "n", main = paste0(gse, " - ", what, " normalization"),
            ylab = "Expression (log2)", cex.main = 1.0, cex.lab = 0.95,
            ylim = range(em, na.rm = TRUE) * c(1, 1.02))
    # group separator + group labels on the x axis (no per-sample labels -> no overlap)
    abline(v = n_ctl + 0.5, lty = 2, col = "grey45")
    mtext(sprintf("Control (n=%d)", n_ctl), side = 1, at = (1 + n_ctl) / 2, line = 1.2, cex = 0.85)
    mtext(sprintf("%s (n=%d)", lbl_case, n_case), side = 1, at = n_ctl + (ncol(em) - n_ctl) / 2 + 0.5, line = 1.2, cex = 0.85)
    axis(2)
    box()
  }
}

svg(file.path(figdir, "Fig1_normalization_boxplots_V1.svg"), width = 14, height = 11)
par(mfrow = c(3, 2), mar = c(4.5, 4.5, 2.6, 1), oma = c(2.5, 0, 2, 0))
for (gse in names(boxdata)) draw_pair(gse, boxdata[[gse]])
legend("topright", inset = c(-0.005, -0.012), xpd = NA, legend = c("Case (DN/DKD)", "Control"),
       fill = c("#E69F00", "#56B4E9"), horiz = TRUE, cex = 0.85, bty = "n")
mtext("Expression distributions before and after quantile normalization (per-sample boxes, grouped)",
      side = 3, outer = TRUE, line = 0.6, cex = 1.0)
dev.off()
cat("Fig1 V1 SVG saved\n")


############################################################
# FILE: step3a_parse_sm.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 3a: Parse GEO series_matrix files -> expression matrix TSV + phenotype CSV.
Outputs per dataset:
  data/{GSE}_expr.tsv   : probe x sample matrix (tab-separated, first col ID_REF)
  data/{GSE}_pheno.csv  : gsm,title,group columns
"""
import gzip, csv, os

BASE = "PROJECT_ROOT/data"
DATASETS = {
    "GSE96804":  dict(in_f="GSE96804_sm.txt.gz",  group_rule="source"),   # DN vs Control by source/title
    "GSE30528":  dict(in_f="GSE30528_sm.txt.gz",  group_rule="disease"),  # DKD vs control
    "GSE30529":  dict(in_f="GSE30529_sm.txt.gz",  group_rule="disease"),
}

def parse(gse, cfg):
    path = os.path.join(BASE, "series_matrix", cfg["in_f"])
    meta = {"gsm": [], "title": [], "source": [], "chars": []}
    expr_lines = []
    in_table = False
    with gzip.open(path, "rt") as f:
        for ln in f:
            if not in_table:
                if ln.startswith("!Sample_geo_accession"):
                    meta["gsm"] = ln.rstrip("\n").split("\t")[1:]
                elif ln.startswith("!Sample_title"):
                    meta["title"] = ln.rstrip("\n").split("\t")[1:]
                elif ln.startswith("!Sample_source_name_ch1"):
                    meta["source"] = ln.rstrip("\n").split("\t")[1:]
                elif ln.startswith("!Sample_characteristics_ch1"):
                    meta["chars"].append(ln.rstrip("\n").split("\t")[1:])
                elif ln.startswith("!series_matrix_table_begin"):
                    in_table = True
                continue
            if ln.startswith("!series_matrix_table_end"):
                break
            expr_lines.append(ln.rstrip("\n"))
    n = len(meta["gsm"])
    assert all(len(meta[k]) == n for k in ["gsm", "title", "source"]) and n > 0, f"meta mismatch {gse}"
    # group assignment
    grp = []
    for i in range(n):
        if cfg["group_rule"] == "source":
            s = meta["source"][i].strip().strip('"')
            t = meta["title"][i].strip().strip('"')
            if "tumor nephrectomies" in s or "Control" in t:
                grp.append("Control")
            else:
                grp.append("DN")
        else:
            d = [c[i].strip().strip('"') for c in meta["chars"] if any("disease state" in c[j] for j in range(len(c)))]
            ds = d[0] if d else ""
            grp.append("DKD" if "DKD" in ds or "diabetic" in ds else "Control")
    # write pheno
    with open(os.path.join(BASE, f"{gse}_pheno.csv"), "w", newline="") as f:
        w = csv.writer(f)
        w.writerow(["gsm", "title", "source", "group"])
        for i in range(n):
            w.writerow([meta["gsm"][i].strip('"'), meta["title"][i].strip('"'),
                        meta["source"][i].strip('"'), grp[i]])
    # write expression
    with open(os.path.join(BASE, f"{gse}_expr.tsv"), "w") as f:
        f.write("\n".join(expr_lines) + "\n")
    from collections import Counter
    print(gse, "samples:", n, "groups:", dict(Counter(grp)),
          "expr rows:", len(expr_lines) - 1, flush=True)

for gse, cfg in DATASETS.items():
    parse(gse, cfg)
print("step3a done")


############################################################
# FILE: step3b_preprocess.R
############################################################
# Step 3b: Preprocess 3 datasets (annotation mapping, probe merging, log2, quantile norm, boxplots)
# Run with: /public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step3b_preprocess.R
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(limma) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step3")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

# ---- platform annotation: probe -> gene symbol ----
ann571 <- read.delim(file.path(base, "data", "GPL571_table.txt"), skip = 1,
                     sep = "\t", quote = "\"", check.names = FALSE, stringsAsFactors = FALSE)
ann571 <- ann571[ann571[[1]] != "!platform_table_end", ]
g571 <- data.frame(probe = ann571$ID,
                   symbol = sub("\\s*$", "", sub("^\\s*", "", sapply(strsplit(as.character(ann571[["Gene Symbol"]]), "///"), function(x) trimws(x[1])))),
                   stringsAsFactors = FALSE)

ann17586 <- read.delim(file.path(base, "data", "GPL17586_table.txt"), skip = 1,
                       sep = "\t", quote = "\"", check.names = FALSE, stringsAsFactors = FALSE)
ann17586 <- ann17586[ann17586[[1]] != "!platform_table_end", ]
ga <- as.character(ann17586[["gene_assignment"]])
first_entry <- sub("\\s*///.*$", "", ga)
sym17586 <- sub("^[^/]*//\\s*", "", first_entry)
sym17586 <- trimws(sub("//.*$", "", paste0(sym17586, "//")))
g17586 <- data.frame(probe = ann17586$ID, symbol = trimws(sym17586), stringsAsFactors = FALSE)

cat("GPL571 probes:", nrow(g571), " with symbol:", sum(g571$symbol != "" & !is.na(g571$symbol)), "\n")
cat("GPL17586 probes:", nrow(g17586), " with symbol:", sum(g17586$symbol != "" & !is.na(g17586$symbol)), "\n")

datasets <- list(
  GSE96804 = list(platform = g17586),
  GSE30528 = list(platform = g571),
  GSE30529 = list(platform = g571)
)

boxdata <- list()
for (gse in names(datasets)) {
  cat("\n====", gse, "====\n")
  expr <- read.delim(file.path(base, "data", paste0(gse, "_expr.tsv")),
                     sep = "\t", header = TRUE, check.names = FALSE, stringsAsFactors = FALSE)
  rownames(expr) <- expr[[1]]; expr[[1]] <- NULL
  expr <- as.matrix(expr)
  storage.mode(expr) <- "numeric"
  pheno <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  stopifnot(identical(as.character(pheno$gsm), colnames(expr)))
  cat("matrix:", nrow(expr), "probes x", ncol(expr), "samples\n")
  # probe -> symbol
  ann <- datasets[[gse]]$platform
  m <- merge(data.frame(probe = rownames(expr), stringsAsFactors = FALSE), ann, by = "probe", sort = FALSE)
  idx <- match(m$probe, rownames(expr))
  em <- expr[idx, , drop = FALSE]
  m$symbol <- ifelse(is.na(m$symbol) | m$symbol == "", paste0(m$probe, "_NA"), m$symbol)
  cat("mapped:", sum(!grepl("_NA$", m$symbol)), "/", nrow(m), "\n")
  rownames(em) <- m$symbol
  em <- avereps(em)   # average duplicated probes per gene
  # log2 decision (mirroring the paper: assess distribution, transform if needed)
  mx <- max(em, na.rm = TRUE); mn <- min(em, na.rm = TRUE)
  if (mx > 100) {
    em <- log2(em + 1); transformed <- TRUE
  } else transformed <- FALSE
  cat(sprintf("range before transform: %.2f - %.2f ; log2 transform applied: %s\n", mn, mx, transformed))
  # quantile normalization (normalizeBetweenArrays default for matrix = quantile)
  em_norm <- normalizeBetweenArrays(em)
  boxdata[[gse]] <- list(before = em, after = em_norm, pheno = pheno)
  out <- data.frame(gene = rownames(em_norm), em_norm, check.names = FALSE)
  write.csv(out, file.path(outdir, paste0(gse, "_norm_expr.csv")), row.names = FALSE)
  write.csv(out, file.path(base, "data", paste0(gse, "_norm_expr.csv")), row.names = FALSE)
  cat("norm matrix saved:", nrow(em_norm), "genes x", ncol(em_norm), "samples\n")
  cat("median before:", round(median(apply(em, 2, median)), 3),
      " after:", round(median(apply(em_norm, 2, median)), 3), "\n")
}

# ---- boxplots (paper Figure 2 style): before vs after for each dataset ----
svg(file.path(outdir, "fig2_normalization_boxplots.svg"), width = 11, height = 10)
par(mfrow = c(3, 2), mar = c(5, 4.5, 3, 1), cex.axis = 0.7)
for (gse in names(boxdata)) {
  em <- boxdata[[gse]]$before; ph <- boxdata[[gse]]$pheno
  ord <- order(ph$group); em <- em[, ord]; lab <- ph$group[ord]
  cols <- ifelse(lab == "DN", "#E69F00", ifelse(lab == "DKD", "#E69F00", "#56B4E9"))
  boxplot(em, outline = FALSE, col = cols, las = 2, pch = ".",
          main = paste0(gse, " before normalization"),
          ylab = "Expression", names = sub("Human glomeruli ", "", ph$title[ord]), cex.main = 0.9)
  boxplot(boxdata[[gse]]$after[, ord], outline = FALSE, col = cols, las = 2, pch = ".",
          main = paste0(gse, " after normalization"),
          ylab = "Expression", names = sub("Human glomeruli ", "", ph$title[ord]), cex.main = 0.9)
}
legend("bottomright", inset = c(-0.02, -0.25), legend = c("Case (DN/DKD)", "Control"),
       fill = c("#E69F00", "#56B4E9"), xpd = TRUE, horiz = TRUE, cex = 0.8)
dev.off()
cat("\nfig2_normalization_boxplots.svg written\n")
cat("step3b done\n")


############################################################
# FILE: step3b_run.sh
############################################################
#!/usr/bin/env bash
# Run step3b preprocessing with the platform R 4.5.3 wrapper
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step3b_preprocess.R


############################################################
# FILE: step4_deg.R
############################################################
# Step 4: DEG analysis in discovery set (GSE96804, 41 DN vs 20 Control) + GRDEGs + Figure 3 panels
# Run with step4_run.sh
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(limma); library(ggplot2); library(pheatmap); library(ggVennDiagram)
})
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step4")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

# ---- load normalized discovery matrix ----
expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL
expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
stopifnot(identical(pheno$gsm, colnames(expr)))
group <- factor(pheno$group, levels = c("Control", "DN"))
cat("dim:", dim(expr), " table:", table(group), "\n")

# ---- limma differential expression (paper: |logFC|>1, BH adj.p<0.05) ----
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
fit <- lmFit(expr, design)
contrast <- makeContrasts(DN - Control, levels = design)
fit2 <- eBayes(contrasts.fit(fit, contrast))
deg <- topTable(fit2, number = Inf, adjust.method = "BH")
deg$gene <- rownames(deg)
deg <- deg[, c("gene", "logFC", "P.Value", "adj.P.Val", "t")]
write.csv(deg, file.path(outdir, "DEG_all_GSE96804.csv"), row.names = FALSE)
sig <- subset(deg, abs(logFC) > 1 & adj.P.Val < 0.05)
sig <- sig[order(sig$adj.P.Val), ]
n_up <- sum(sig$logFC > 1); n_down <- sum(sig$logFC < -1)
cat("DEGs:", nrow(sig), " up:", n_up, " down:", n_down, "\n")
sig$direction <- ifelse(sig$logFC > 0, "up", "down")
write.csv(sig, file.path(outdir, "DEG_sig_GSE96804.csv"), row.names = FALSE)

# ---- GRDEGs: intersect with GRG list ----
grg <- read.csv(file.path(base, "output", "GRG_list.csv"), stringsAsFactors = FALSE)
cat("GRG list:", nrow(grg), "\n")
grdegs <- sig[sig$gene %in% grg$gene_symbol, ]
grdegs <- grdegs[order(grdegs$adj.P.Val), ]
cat("GRDEGs:", nrow(grdegs), " up:", sum(grdegs$logFC > 0), " down:", sum(grdegs$logFC < 0), "\n")
write.csv(grdegs, file.path(outdir, "GRDEGs_GSE96804.csv"), row.names = FALSE)
sig_genes <- grdegs$gene
writeLines(sig_genes, file.path(outdir, "GRDEG_symbols.txt"))

# ---- Fig 3A volcano ----
ggrepel_free <- function(df) {
  geom_text(data = df, aes(x = logFC, y = -log10(adj.P.Val), label = gene),
            size = 2.3, vjust = -0.7, color = "black")
}
print_plot <- function(p, path, w, h) {
  ggsave(path, p, width = w, height = h, device = grDevices::svg)
  cat("saved:", path, "\n")
}
p <- ggplot(as.data.frame(deg), aes(x = logFC, y = -log10(adj.P.Val))) +
  geom_point(aes(color = ifelse(gene %in% sig_genes, "GRDEG",
                          ifelse(abs(logFC) > 1 & adj.P.Val < 0.05, "DEG", "NS"))),
             size = 1.1, alpha = 0.6) +
  scale_color_manual(values = c("GRDEG" = "#D62728", "DEG" = "#7F7F7F", "NS" = "#CCCCCC"),
                     name = NULL) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey40") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey40") +
  labs(title = "Volcano plot: DN vs Control (GSE96804)",
       subtitle = paste0(nrow(sig), " DEGs (|logFC|>1, adj.p<0.05); ",
                         nrow(grdegs), " GRDEGs highlighted"),
       x = "log2 Fold Change", y = "-log10(adjusted p-value)") +
  theme_bw(base_size = 11) + theme(legend.position = "right")
top_lab <- head(sig[order(-abs(sig$logFC)), ], 12)
p <- p + ggrepel_free(top_lab)
print_plot(p, file.path(outdir, "fig3a_volcano.svg"), 9, 7)

# ---- Fig 3B Venn ----
svg(file.path(outdir, "fig3b_venn.svg"), width = 7, height = 7)
ggVennDiagram(list("DEGs (GSE96804)" = sig$gene, "GRGs" = grg$gene),
              label = "count", edge_color = "grey30") +
  scale_fill_gradient(low = "#F4FAFF", high = "#D62728") +
  ggtitle(paste0("DEGs n GRGs = ", nrow(grdegs), " GRDEGs"))
dev.off(); cat("saved fig3b_venn.svg\n")

# ---- Fig 3C heatmap top 20 GRDEGs ----
top20 <- head(grdegs, 20)$gene
top20 <- intersect(top20, rownames(expr))
mat <- expr[top20, , drop = FALSE]
mat <- t(scale(t(mat)))
ann_col <- data.frame(Group = group); rownames(ann_col) <- colnames(expr)
svg(file.path(outdir, "fig3c_heatmap_top20.svg"), width = 12, height = 7)
pheatmap(mat, scale = "none", annotation_col = ann_col,
         annotation_colors = list(Group = c(Control = "#56B4E9", DN = "#E69F00")),
         color = colorRampPalette(c("#3B4CC0", "white", "#B40426"))(100),
         cluster_cols = TRUE, show_colnames = FALSE,
         main = "Top 20 GRDEGs (GSE96804)", fontsize_row = 8)
dev.off(); cat("saved fig3c_heatmap_top20.svg\n")
writeLines(top20, file.path(outdir, "top20_GRDEG_symbols.txt"))

# ---- Fig 3D chromosomal localization (RCircos; fallback bar) ----
ok_circos <- tryCatch({
  suppressPackageStartupMessages({ library(RCircos); library(org.Hs.eg.db); library(AnnotationDbi) })
  sym2eg <- toTable(org.Hs.egALIAS2EG)
  loc <- toTable(org.Hs.egCHRLOC); chr <- toTable(org.Hs.egCHR)
  ids <- sym2eg[sym2eg$alias_symbol %in% sig_genes, ]
  ids <- ids[!duplicated(ids$alias_symbol), ]
  m <- merge(ids, loc, by.x = "gene_id", by.y = "gene_id")
  m <- merge(m, chr[, c("gene_id", "chromosome")], by = "gene_id")
  m <- m[m$chromosome %in% as.character(1:22), ]
  m$loc <- abs(m$start_position)
  m <- m[!duplicated(m$alias_symbol), ]
  m$value <- abs(grdegs$logFC[match(m$alias_symbol, grdegs$gene)])
  plot_data <- data.frame(Chromosome = paste0("chr", m$chromosome),
                          chromStart = m$loc, chromEnd = m$loc + 1e5,
                          value = m$value,
                          gene = m$alias_symbol, stringsAsFactors = FALSE)
  plot_data <- plot_data[order(plot_data$Chromosome, plot_data$chromStart), ]
  data(UCSC.HG38.Human.CytoBandIdeogram, package = "RCircos")
  cyto <- UCSC.HG38.Human.CytoBandIdeogram
  cyto <- cyto[!cyto$Chromosome %in% c("chrX", "chrY"), ]
  svg(file.path(outdir, "fig3d_circos.svg"), width = 9, height = 9)
  RCircos.Set.Core.Components(cyto, chr.exclude = NULL, tracks.inside = 4, tracks.outside = 0)
  RCircos.Chromosome.Ideogram.Plot()
  RCircos.Scatter.Plot(plot_data, data.col = 4, track.num = 1, side = "in", by.fold = 10)
  RCircos.Gene.Name.Plot(plot_data, name.col = 5, track.num = 2, side = "in")
  dev.off()
  cat("saved fig3d_circos.svg (", nrow(plot_data), "GRDEGs mapped )\n")
  TRUE
}, error = function(e) { cat("RCircos failed:", conditionMessage(e), "\n"); FALSE })

if (!ok_circos || !file.exists(file.path(outdir, "fig3d_circos.svg"))) {
  # fallback: chromosome distribution bar plot
  suppressPackageStartupMessages({ library(org.Hs.eg.db); library(AnnotationDbi) })
  sym2eg <- toTable(org.Hs.egALIAS2EG); chr <- toTable(org.Hs.egCHR)
  ids <- sym2eg[sym2eg$alias_symbol %in% sig_genes, ]
  m <- merge(ids, chr[, c("gene_id", "chromosome")], by = "gene_id")
  cnt <- sort(table(m$chromosome), decreasing = TRUE)
  df <- data.frame(chr = names(cnt), n = as.integer(cnt))
  p <- ggplot(df, aes(x = reorder(chr, -n), y = n)) +
    geom_col(fill = "#D62728") + theme_bw() +
    labs(title = "GRDEGs chromosomal distribution (GSE96804)",
         x = "Chromosome", y = "Number of GRDEGs")
  ggsave(file.path(outdir, "fig3d_chr_distribution.svg"), p, width = 8, height = 4.5,
         device = grDevices::svg)
  cat("saved fallback fig3d_chr_distribution.svg\n")
}

cat("\nSTEP4 SUMMARY: DEGs:", nrow(sig), "(up", n_up, "/ down", n_down,
    ")  GRDEGs:", nrow(grdegs), "\nTop20:", paste(top20, collapse = ", "), "\n")


############################################################
# FILE: step4_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step4_deg.R


############################################################
# FILE: step4d_circos.R
############################################################
# Step 4d: RCircos chromosomal localization of GRDEGs (retry with NA removal)
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({ library(RCircos); library(org.Hs.eg.db); library(AnnotationDbi) })
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step4")
grdegs <- read.csv(file.path(outdir, "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
sym2eg <- toTable(org.Hs.egALIAS2EG)
loc <- toTable(org.Hs.egCHRLOC)   # gene_id, start_location, Chromosome
ids <- sym2eg[sym2eg$alias_symbol %in% grdegs$gene, ]
ids <- ids[!duplicated(ids$alias_symbol), ]
m <- merge(ids, loc, by = "gene_id")
m <- m[m$Chromosome %in% as.character(1:22), ]
m$loc <- abs(as.numeric(as.character(m$start_location)))
m <- m[!is.na(m$loc) & m$loc > 0, ]
m <- m[!duplicated(m$alias_symbol), ]
m$value <- abs(grdegs$logFC[match(m$alias_symbol, grdegs$gene)])
m <- m[!is.na(m$value), ]
plot_data <- data.frame(Chromosome = paste0("chr", m$Chromosome),
                        chromStart = m$loc, chromEnd = m$loc + 1e5,
                        value = m$value, gene = m$alias_symbol, stringsAsFactors = FALSE)
plot_data <- plot_data[order(plot_data$Chromosome, plot_data$chromStart), ]
cat("GRDEGs mapped to chromosomes:", nrow(plot_data), "\n")
data(UCSC.HG38.Human.CytoBandIdeogram, package = "RCircos")
cyto <- UCSC.HG38.Human.CytoBandIdeogram
cyto <- cyto[!cyto$Chromosome %in% c("chrX", "chrY"), ]
svg(file.path(outdir, "fig3d_circos.svg"), width = 9, height = 9)
RCircos.Set.Core.Components(cyto, chr.exclude = NULL, tracks.inside = 4, tracks.outside = 0)
plot.new()
RCircos.Chromosome.Ideogram.Plot()
RCircos.Scatter.Plot(plot_data, data.col = 4, track.num = 1, side = "in", by.fold = 10)
RCircos.Gene.Name.Plot(plot_data, name.col = 5, track.num = 2, side = "in")
dev.off()
cat("saved fig3d_circos.svg\n")
# chromosome count table for the report
cnt <- table(plot_data$Chromosome)
cnt <- cnt[order(as.integer(sub("chr", "", names(cnt))), decreasing = FALSE)]
write.table(data.frame(chromosome = names(cnt), n_grdegs = as.integer(cnt)),
            file.path(outdir, "GRDEGs_chromosome_counts.csv"),
            sep = ",", row.names = FALSE, quote = FALSE)
print(cnt)


############################################################
# FILE: step4d_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step4d_circos.R


############################################################
# FILE: step5_enrich.R
############################################################
# Step 5: GO/KEGG enrichment of 28 GRDEGs (paper Fig 4: bubble + network plots)
# Thresholds mirroring the paper: adj.p<0.05 (BH) and qvalue<0.25
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(clusterProfiler); library(org.Hs.eg.db); library(enrichplot)
  library(ggplot2); library(gridExtra)
})
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step5")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

grdegs <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
genes <- grdegs$gene
fc <- grdegs$logFC; names(fc) <- genes
cat("GRDEGs:", length(genes), "\n")

eg <- bitr(genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
cat("mapped to ENTREZID:", nrow(eg), "/", length(genes), "\n")
gene_ids <- eg$ENTREZID
fc_id <- fc[match(eg$SYMBOL, genes)]; names(fc_id) <- eg$ENTREZID

# ---- GO enrichment (local org.Hs.eg.db) ----
res <- list()
for (ont in c("BP", "CC", "MF")) {
  e <- enrichGO(gene = gene_ids, OrgDb = org.Hs.eg.db, keyType = "ENTREZID",
                ont = ont, pvalueCutoff = 0.05, qvalueCutoff = 0.25,
                pAdjustMethod = "BH", readable = TRUE)
  cat("GO", ont, "significant terms:", if (is.null(e)) 0 else nrow(as.data.frame(e)), "\n")
  res[[ont]] <- e
  if (!is.null(e) && nrow(as.data.frame(e)) > 0)
    write.csv(as.data.frame(e), file.path(outdir, paste0("GO_", ont, ".csv")), row.names = FALSE)
}

# ---- KEGG enrichment via offline TERM2GENE (from rest.kegg.jp dump) ----
link <- read.delim(file.path(base, "data", "kegg_hsa_link.tsv"), header = FALSE,
                   col.names = c("path", "gene"), stringsAsFactors = FALSE)
link$pathway <- sub("path:", "", link$path)
link$entrez <- sub("hsa:", "", link$gene)
pw <- read.delim(file.path(base, "data", "kegg_hsa_pathways.tsv"), header = FALSE,
                 stringsAsFactors = FALSE)
names(pw) <- c("pathway", "name")
pw$name <- sub(" - Homo sapiens.*$", "", pw$name)
t2g <- unique(link[, c("pathway", "entrez")])
t2n <- pw[, c("pathway", "name")]
kegg <- enricher(gene = gene_ids, TERM2GENE = t2g, TERM2NAME = t2n,
                 pvalueCutoff = 0.05, qvalueCutoff = 0.25, pAdjustMethod = "BH")
if (!is.null(kegg) && nrow(as.data.frame(kegg)) > 0) {
  kegg <- setReadable(kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
  write.csv(as.data.frame(kegg), file.path(outdir, "KEGG.csv"), row.names = FALSE)
  cat("KEGG significant terms:", nrow(as.data.frame(kegg)), "\n")
} else {
  cat("KEGG: no significant terms\n")
}
res[["KEGG"]] <- kegg

# ---- Fig 4A style bubble plot: top10 per category in 2x2 grid ----
mk_bubble <- function(e, ttl) {
  if (is.null(e) || nrow(as.data.frame(e)) == 0) {
    return(nullGROB())
  }
  d <- as.data.frame(e)[1:min(10, nrow(as.data.frame(e))), ]
  d$GeneRatio <- sapply(d$GeneRatio, function(r) { a <- as.numeric(strsplit(r, "/")[[1]]); a[1]/a[2] })
  d$Description <- factor(d$Description, levels = rev(d$Description))
  ggplot(d, aes(x = GeneRatio, y = Description)) +
    geom_point(aes(size = Count, color = -log10(p.adjust))) +
    scale_color_gradient(low = "#FEE0D2", high = "#67000D", name = "-log10(p.adj)") +
    scale_size(range = c(2, 6)) +
    labs(title = ttl, x = "Gene Ratio", y = NULL) +
    theme_bw(base_size = 9) + theme(axis.text.y = element_text(size = 7),
                                    plot.title = element_text(size = 10, face = "bold"))
}
nullGROB <- function() grid::nullGrob()
plots <- list(
  mk_bubble(res$BP, "GO Biological Process"),
  mk_bubble(res$CC, "GO Cellular Component"),
  mk_bubble(res$MF, "GO Molecular Function"),
  mk_bubble(res$KEGG, "KEGG Pathway")
)
g <- gridExtra::arrangeGrob(grobs = plots, ncol = 2, nrow = 2)
ggsave(file.path(outdir, "fig4a_enrichment_bubble.svg"), g, width = 14, height = 9,
       device = grDevices::svg)
cat("saved fig4a_enrichment_bubble.svg\n")

# ---- Fig 4B-E style network plots (cnetplot, top5 terms each) ----
for (nm in names(res)) {
  e <- res[[nm]]
  if (is.null(e) || nrow(as.data.frame(e)) == 0) next
  svg(file.path(outdir, paste0("fig4_cnetplot_", nm, ".svg")), width = 9, height = 8)
  tryCatch({
    print(cnetplot(e, showCategory = 5, foldChange = fc_id, cex.params = list(gene_label = 0.55, category_label = 0.8)) +
            ggplot2::ggtitle(paste0("GRDEGs ", nm, " network (top 5 terms)")) +
            ggplot2::theme(plot.title = ggplot2::element_text(size = 11)))
  }, error = function(err) {
    cat("cnetplot", nm, "failed:", conditionMessage(err), "\n")
    plot.new()
  })
  dev.off()
  cat("saved fig4_cnetplot_", nm, ".svg\n", sep = "")
}

# summary
cat("\nSTEP5 SUMMARY\n")
for (nm in names(res)) {
  e <- res[[nm]]
  cat(nm, ":", if (is.null(e) || nrow(as.data.frame(e)) == 0) 0 else nrow(as.data.frame(e)), "terms\n")
}


############################################################
# FILE: step5_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step5_enrich.R


############################################################
# FILE: step5b_cnetplot.R
############################################################
# Step 5b: cnetplot retry with correct enrichplot parameter names
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(clusterProfiler); library(org.Hs.eg.db); library(enrichplot); library(ggplot2)
})
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step5")
grdegs <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
fc <- grdegs$logFC; names(fc) <- grdegs$gene
eg <- bitr(grdegs$gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
fc_id <- fc[match(eg$SYMBOL, grdegs$gene)]; names(fc_id) <- eg$ENTREZID

redo_go <- function(ont) enrichGO(gene = eg$ENTREZID, OrgDb = org.Hs.eg.db, keyType = "ENTREZID",
                                  ont = ont, pvalueCutoff = 0.05, qvalueCutoff = 0.25,
                                  pAdjustMethod = "BH", readable = TRUE)
link <- read.delim(file.path(base, "data", "kegg_hsa_link.tsv"), header = FALSE,
                   col.names = c("path", "gene"))
link$pathway <- sub("path:", "", link$path); link$entrez <- sub("hsa:", "", link$gene)
pw <- read.delim(file.path(base, "data", "kegg_hsa_pathways.tsv"), header = FALSE)
names(pw) <- c("pathway", "name"); pw$name <- sub(" - Homo sapiens.*$", "", pw$name)
kegg <- enricher(gene = eg$ENTREZID, TERM2GENE = unique(link[, c("pathway", "entrez")]),
                 TERM2NAME = pw[, c("pathway", "name")], pvalueCutoff = 0.05, qvalueCutoff = 0.25)
if (!is.null(kegg) && nrow(as.data.frame(kegg)) > 0)
  kegg <- setReadable(kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")

objs <- list(BP = redo_go("BP"), CC = redo_go("CC"), MF = redo_go("MF"), KEGG = kegg)
for (nm in names(objs)) {
  e <- objs[[nm]]
  if (is.null(e) || nrow(as.data.frame(e)) == 0) next
  svg(file.path(outdir, paste0("fig4_cnetplot_", nm, ".svg")), width = 12, height = 10)
  p <- try(cnetplot(e, showCategory = 5, foldChange = fc_id,
                    size_category = 1, node_label = "all"), silent = TRUE)
  if (inherits(p, "try-error")) {
    cat("cnetplot", nm, "still failed:", as.character(p), "\n")
    plot.new(); text(0.5, 0.5, paste(nm, "cnetplot failed"))
  } else {
    print(p + ggtitle(paste0("GRDEGs ", nm, " network (top 5 terms)")) +
            theme(plot.title = element_text(size = 11)))
  }
  dev.off()
  cat("saved fig4_cnetplot_", nm, ".svg\n", sep = "")
}


############################################################
# FILE: step5b_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step5b_cnetplot.R


############################################################
# FILE: step6_model.R
############################################################
# Step 6: Molecular recognition model (univariate logistic -> RF -> LASSO), paper Fig 5
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(glmnet); library(randomForest); library(ggplot2)
})
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step6")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
stopifnot(identical(pheno$gsm, colnames(expr)))
grdegs <- read.csv(file.path(base, "output", "step4", "GRDEGs_GSE96804.csv"), stringsAsFactors = FALSE)
gset <- intersect(grdegs$gene, rownames(expr))
X <- t(expr[gset, ])                    # samples x genes
y <- factor(pheno$group, levels = c("Control", "DN"))
cat("X:", dim(X), " y:", table(y), "\n")

# ---- Step 1: univariate logistic regression (p<0.05), forest plot ----
uni <- do.call(rbind, lapply(gset, function(g) {
  d <- data.frame(y = y, x = X[, g])
  fit <- suppressWarnings(glm(y ~ x, family = binomial, data = d))
  s <- summary(fit)$coefficients
  est <- s["x", "Estimate"]; se <- s["x", "Std. Error"]; p <- s["x", "Pr(>|z|)"]
  data.frame(gene = g, OR = exp(est), CI_low = exp(est - 1.96 * se),
             CI_high = exp(est + 1.96 * se), coef = est, p = p, stringsAsFactors = FALSE)
}))
uni <- uni[order(uni$p), ]
uni_sig <- subset(uni, p < 0.05)
cat("univariate logistic significant genes:", nrow(uni_sig), "\n")
write.csv(uni, file.path(outdir, "univariate_logistic.csv"), row.names = FALSE)
write.csv(uni_sig, file.path(outdir, "univariate_logistic_sig.csv"), row.names = FALSE)

# Fig 5A forest plot
d <- uni_sig
d$gene <- factor(d$gene, levels = rev(d$gene))
pA <- ggplot(d, aes(x = OR, y = gene)) +
  geom_point(size = 2.2, color = "#2166AC") +
  geom_errorbarh(aes(xmin = pmax(CI_low, 1e-3), xmax = CI_high), height = 0.22, color = "#2166AC") +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey50") +
  scale_x_log10() +
  labs(title = "Univariate logistic regression (GRDEGs, p<0.05)",
       subtitle = paste0(nrow(d), " significant genes; OR on log scale (DN vs Control)"),
       x = "Odds Ratio (log10)", y = NULL) +
  theme_bw(base_size = 10)
ggsave(file.path(outdir, "fig5a_forest_univariate.svg"), pA, width = 8, height = 5.5,
       device = grDevices::svg)
cat("saved fig5a_forest_univariate.svg\n")

# ---- Step 2: random forest (ntree=1000, seed 2024), MeanDecreaseGini top10 ----
set.seed(2024)
rf <- randomForest(x = X, y = y, ntree = 1000, importance = TRUE)
imp <- importance(rf, type = 2, scale = FALSE)   # MeanDecreaseGini
imp <- data.frame(gene = rownames(imp), MeanDecreaseGini = as.numeric(imp[, "MeanDecreaseGini"]))
imp <- imp[order(-imp$MeanDecreaseGini), ]
cat("RF top10 (MeanDecreaseGini):\n"); print(head(imp, 10))
write.csv(imp, file.path(outdir, "RF_importance_MeanDecreaseGini.csv"), row.names = FALSE)
top10 <- head(imp, 10)$gene
writeLines(top10, file.path(outdir, "RF_top10_genes.txt"))
cat("RF OOB error rate:", tail(rf$err.rate[, "OOB"], 1), "\n")

# Fig 5B MeanDecreaseGini scatter (all GRDEGs, top10 highlighted)
imp$top10 <- imp$gene %in% top10
pB <- ggplot(imp, aes(x = reorder(gene, MeanDecreaseGini), y = MeanDecreaseGini, color = top10)) +
  geom_point(size = 2.4) +
  scale_color_manual(values = c("FALSE" = "grey65", "TRUE" = "#B2182B"), guide = "none") +
  labs(title = "Random forest feature importance (MeanDecreaseGini)",
       subtitle = "Top 10 genes in red; GRDEGs ordered ascending",
       x = NULL, y = "MeanDecreaseGini") +
  coord_flip() + theme_bw(base_size = 10) +
  theme(axis.text.y = element_text(size = 6.5))
ggsave(file.path(outdir, "fig5b_rf_importance.svg"), pB, width = 7.5, height = 7,
       device = grDevices::svg)
cat("saved fig5b_rf_importance.svg\n")

# ---- Step 3: LASSO on top10 (cv.glmnet), final model genes ----
X10 <- X[, top10]
cvfit <- cv.glmnet(x = X10, y = y, family = "binomial", alpha = 1, type.measure = "auc",
                   nfolds = 10, standardize = TRUE)
lam <- cvfit$lambda.min
coef <- as.matrix(coef(cvfit, s = "lambda.min"))
nz <- data.frame(gene = rownames(coef), coefficient = as.numeric(coef))
nz <- nz[nz$coefficient != 0 & nz$gene != "(Intercept)", ]
cat("LASSO final model genes (lambda.min =", signif(lam, 4), "):\n"); print(nz)
write.csv(nz, file.path(outdir, "LASSO_model_genes.csv"), row.names = FALSE)
model_genes <- nz$gene
writeLines(model_genes, file.path(outdir, "model_genes.txt"))
write.csv(data.frame(model_genes = model_genes, coefficient = nz$coefficient),
          file.path(outdir, "model_risk_score_formula.csv"), row.names = FALSE)
formula_txt <- paste0("RiskScore = ", paste0(
  sprintf("(%.6f x Exp(%s))", nz$coefficient, model_genes), collapse = " + "))
writeLines(formula_txt, file.path(outdir, "risk_score_formula.txt"))
cat(formula_txt, "\n")

# Fig 5C: lambda selection (cv curve) & Fig 5D: coefficient trajectory
svg(file.path(outdir, "fig5c_lasso_cv.svg"), width = 9.5, height = 5.5)
par(mar = c(4.5, 4.5, 3, 1))
plot(cvfit)
title("LASSO 10-fold cross-validation (binomial AUC)", line = 2.2, cex.main = 0.95)
abline(v = log(lam), lty = 2, col = "#B2182B")
dev.off(); cat("saved fig5c_lasso_cv.svg\n")

svg(file.path(outdir, "fig5d_lasso_trajectory.svg"), width = 9.5, height = 5.5)
par(mar = c(4.5, 4.5, 3, 7))
fitfull <- glmnet(x = X10, y = y, family = "binomial", alpha = 1, standardize = TRUE)
plot(fitfull, xvar = "lambda", label = FALSE)
cf10 <- as.matrix(coef(fitfull, s = lam))
cf10 <- cf10[rownames(cf10) != "(Intercept)", ]
ord10 <- order(cf10); ys10 <- as.numeric(cf10[ord10]); lab10 <- rownames(cf10)[ord10]
minsep <- max(diff(range(ys10)) * 0.055, 1e-3)
for (k in 2:length(ys10)) if (ys10[k] - ys10[k-1] < minsep) ys10[k] <- ys10[k-1] + minsep
lg10 <- log(fitfull$lambda); xr10 <- diff(range(lg10))
text(x = max(lg10) + xr10 * 0.02, y = ys10, labels = lab10, cex = 0.62, adj = 0, xpd = NA)
title("LASSO coefficient trajectories (top10 GRDEGs)", line = 2.2, cex.main = 0.95)
abline(v = log(lam), lty = 2, col = "#B2182B")
dev.off(); cat("saved fig5d_lasso_trajectory.svg\n")

cat("\nSTEP6 SUMMARY\n")
cat("univariate sig:", nrow(uni_sig), " RF top10:", paste(top10, collapse = ","),
    "\nfinal model genes:", paste(model_genes, collapse = ","), "\n")


############################################################
# FILE: step6_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step6_model.R


############################################################
# FILE: step7_evaluate.R
############################################################
# Step 7: Training-set evaluation of the 7-gene DN model (paper Fig 6 A-E)
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(rms); library(pROC); library(ggplot2); library(pheatmap)
  library(GOSemSim); library(org.Hs.eg.db); library(clusterProfiler)
})
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step7")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

expr <- read.csv(file.path(base, "data", "GSE96804_norm_expr.csv"),
                 check.names = FALSE, stringsAsFactors = FALSE)
rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
pheno <- read.csv(file.path(base, "data", "GSE96804_pheno.csv"), stringsAsFactors = FALSE)
stopifnot(identical(pheno$gsm, colnames(expr)))
mg <- readLines(file.path(base, "output", "step6", "model_genes.txt"))
coefs <- read.csv(file.path(base, "output", "step6", "LASSO_model_genes.csv"), stringsAsFactors = FALSE)
cf <- coefs$coefficient; names(cf) <- coefs$gene
stopifnot(setequal(mg, coefs$gene))
y <- factor(pheno$group, levels = c("Control", "DN"))
df <- as.data.frame(t(expr[mg, ]))
df$y <- y
cat("model genes:", paste(mg, collapse = ","), "\n")

# RiskScore (paper formula: sum of coefficient x expression)
rs <- as.numeric(as.matrix(df[, mg]) %*% cf[mg])
df$RiskScore <- rs
write.csv(data.frame(gsm = pheno$gsm, group = pheno$group, RiskScore = rs),
          file.path(outdir, "training_riskscore.csv"), row.names = FALSE)

# ---- Fig 6A nomogram (rms) ----
dd <- datadist(df); options(datadist = "dd")
f <- lrm(y ~ RiskScore, data = df)
svg(file.path(outdir, "fig6a_nomogram.svg"), width = 11, height = 6)
plot(nomogram(f, fun = plogis, funlabel = "P(DN)",
     fun.at = c(0.01, 0.5, 0.99), nint = 5),
     main = "Nomogram: DN molecular recognition model (GSE96804)", cex.axis = 0.75)
dev.off(); cat("saved fig6a_nomogram.svg\n")

# logistic model on 7 genes for predicted probability (for calibration/DCA)
f7 <- lrm(y ~ ., data = df[, c(mg, "y")], x = TRUE, y = TRUE)
pred <- predict(f7, type = "fitted")

# save training intercept so validation nomograms can reuse the frozen model
writeLines(sprintf("%.10f", f7$coefficients[1]),
           file.path(outdir, "model_intercept.txt"))

# ---- Fig 6B ROC / AUC (RiskScore) ----
roc_obj <- roc(response = y, predictor = rs, quiet = TRUE)
auc_val <- auc(roc_obj)
ci <- ci.auc(roc_obj)
cat(sprintf("Training AUC (RiskScore) = %.4f  (95%%CI %.4f-%.4f)\n", auc_val, ci[1], ci[3]))
svg(file.path(outdir, "fig6b_roc_training.svg"), width = 6.5, height = 6)
plot(roc_obj, col = "#B2182B", lwd = 2.2, main = "ROC: DN molecular recognition model (GSE96804)",
     print.auc = TRUE, print.auc.cex = 0.9, cex.main = 0.95, cex.lab = 0.95)
grid()
dev.off(); cat("saved fig6b_roc_training.svg\n")
writeLines(data.frame(AUC = auc_val, CI_low = ci[1], CI_high = ci[3]) |>
             sprintf(fmt = "%.4f"), file.path(outdir, "training_auc.txt"))

# ---- Fig 6C calibration curve (bootstrap) ----
svg(file.path(outdir, "fig6c_calibration.svg"), width = 7, height = 7)
par(mar = c(9.5, 4.5, 3, 2))
cal <- calibrate(f7, method = "boot", B = 300)
plot(cal, xlab = "Predicted probability of DN", ylab = "Observed proportion (DN)",
     main = "Calibration curve (bootstrap B=300)", cex.main = 0.95, cex.lab = 0.95)
abline(0, 1, lty = 2, col = "grey50")
dev.off(); cat("saved fig6c_calibration.svg\n")

# ---- Fig 6D decision curve analysis (net benefit, self-implemented) ----
thr <- seq(0.01, 0.99, by = 0.01)
yy <- as.integer(y == "DN")
n <- length(yy)
prev <- mean(yy)
nb_model <- sapply(thr, function(t) {
  tp <- sum(pred >= t & yy == 1); fp <- sum(pred >= t & yy == 0)
  (tp - fp * (t / (1 - t))) / n
})
nb_all <- prev - (1 - prev) * (thr / (1 - thr))
dca <- data.frame(threshold = thr,
                  NB = c(nb_model, nb_all, rep(0, length(thr))),
                  Strategy = rep(c("Model", "Treat all", "Treat none"), each = length(thr)))
write.csv(data.frame(threshold = thr, NB_model = nb_model, NB_all = nb_all, NB_none = 0),
          file.path(outdir, "training_DCA_netbenefit.csv"), row.names = FALSE)
pD <- ggplot(dca, aes(x = threshold, y = NB, color = Strategy)) +
  geom_line(linewidth = 0.9) +
  coord_cartesian(ylim = c(-0.05, max(nb_model) * 1.05)) +
  labs(title = "Decision curve analysis (training set GSE96804)",
       x = "Threshold probability", y = "Net benefit") +
  theme_bw(base_size = 10) + theme(legend.position = "right")
ggsave(file.path(outdir, "fig6d_dca_training.svg"), pD, width = 7, height = 5,
       device = grDevices::svg)
cat("saved fig6d_dca_training.svg\n")

# ---- Fig 6E GO semantic similarity (friends analysis, GOSemSim) ----
eg <- bitr(mg, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
cat("genes with ENTREZID:", nrow(eg), "\n")
semBP <- godata(org.Hs.eg.db, ont = "BP")
sim <- try(mgeneSim(genes = eg$ENTREZID, semData = semBP, measure = "Wang"), silent = TRUE)
if (!inherits(sim, "try-error") && !is.null(sim)) {
  sym <- eg$SYMBOL[match(colnames(sim), eg$ENTREZID)]
  colnames(sim) <- sym; rownames(sim) <- sym
  sim <- sim[complete.cases(sim), complete.cases(sim)]
  svg(file.path(outdir, "fig6e_gosemsim_friends.svg"), width = 7, height = 6)
  pheatmap(sim, display_numbers = TRUE, number_format = "%.2f", number_size = 7,
           color = colorRampPalette(c("#FFF5F0", "#67000D"))(100),
           main = "GO(BP) semantic similarity of model genes (friends)")
  dev.off()
  avg <- sort(rowMeans(sim, na.rm = TRUE), decreasing = TRUE)
  write.csv(data.frame(gene = names(avg), mean_similarity = as.numeric(avg)),
            file.path(outdir, "friends_mean_similarity.csv"), row.names = FALSE)
  cat("friends (mean GO-BP similarity) ranking:\n"); print(avg)
} else {
  cat("GOSemSim failed:", as.character(sim), "\n")
}

# RiskScore group comparison (extra support)
p <- ggplot(data.frame(group = y, RiskScore = rs), aes(x = group, y = RiskScore, fill = group)) +
  geom_boxplot(outlier.size = 0.8) + geom_jitter(width = 0.15, size = 0.7, alpha = 0.5) +
  scale_fill_manual(values = c(Control = "#56B4E9", DN = "#E69F00"), guide = "none") +
  labs(title = "RiskScore by group (training set)", x = NULL, y = "RiskScore") +
  theme_bw(base_size = 10)
ggsave(file.path(outdir, "fig6f_riskscore_box.svg"), p, width = 5.5, height = 4.5,
       device = grDevices::svg)
cat("saved fig6f_riskscore_box.svg\n")
wt <- wilcox.test(RiskScore ~ y, data = data.frame(RiskScore = rs, y = y))$p.value
cat("RiskScore group Wilcoxon p =", signif(wt, 4), "\n")
cat("\nSTEP7 SUMMARY: AUC =", sprintf("%.4f", auc_val), "\n")


############################################################
# FILE: step7_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step7_evaluate.R


############################################################
# FILE: step7b_friends.R
############################################################
# Step 7b: GOSemSim friends analysis + RiskScore box (retry with clusterProfiler loaded)
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(clusterProfiler); library(GOSemSim); library(org.Hs.eg.db)
  library(pheatmap); library(ggplot2)
})
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step7")
mg <- readLines(file.path(base, "output", "step6", "model_genes.txt"))
rs <- read.csv(file.path(outdir, "training_riskscore.csv"), stringsAsFactors = FALSE)

eg <- bitr(mg, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
cat("genes with ENTREZID:", nrow(eg), "\n")
semBP <- godata(org.Hs.eg.db, ont = "BP")
sim <- try(mgeneSim(genes = eg$ENTREZID, semData = semBP, measure = "Wang"), silent = TRUE)
if (!inherits(sim, "try-error") && !is.null(sim)) {
  sym <- eg$SYMBOL[match(colnames(sim), eg$ENTREZID)]
  colnames(sim) <- sym; rownames(sim) <- sym
  keep <- complete.cases(sim) & colSums(!is.na(sim)) > 0
  sim <- sim[keep, keep]
  svg(file.path(outdir, "fig6e_gosemsim_friends.svg"), width = 7, height = 6)
  pheatmap(sim, display_numbers = TRUE, number_format = "%.2f", number_size = 7,
           color = colorRampPalette(c("#FFF5F0", "#67000D"))(100),
           main = "GO(BP) semantic similarity of model genes (friends)")
  dev.off()
  avg <- sort(rowMeans(sim, na.rm = TRUE), decreasing = TRUE)
  write.csv(data.frame(gene = names(avg), mean_similarity = as.numeric(avg)),
            file.path(outdir, "friends_mean_similarity.csv"), row.names = FALSE)
  cat("friends (mean GO-BP similarity) ranking:\n"); print(avg)
} else {
  cat("GOSemSim failed:", as.character(sim), "\n")
}

p <- ggplot(rs, aes(x = group, y = RiskScore, fill = group)) +
  geom_boxplot(outlier.size = 0.8) + geom_jitter(width = 0.15, size = 0.7, alpha = 0.5) +
  scale_fill_manual(values = c(Control = "#56B4E9", DN = "#E69F00"), guide = "none") +
  labs(title = "RiskScore by group (training set)", x = NULL, y = "RiskScore") +
  theme_bw(base_size = 10)
ggsave(file.path(outdir, "fig6f_riskscore_box.svg"), p, width = 5.5, height = 4.5,
       device = grDevices::svg)
rs$group <- factor(rs$group, levels = c("Control", "DN"))
cat("RiskScore group Wilcoxon p =", signif(wilcox.test(RiskScore ~ group, data = rs)$p.value, 4), "\n")
cat("saved fig6f_riskscore_box.svg\n")


############################################################
# FILE: step7b_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step7b_friends.R


############################################################
# FILE: step8_run.sh
############################################################
#!/usr/bin/env bash
set -e
cd PROJECT_ROOT/scripts
/public/tools/r-packages/R-4.5.3/R --vanilla --slave -f step8_validate.R


############################################################
# FILE: step8_validate.R
############################################################
# Step 8: External validation (GSE30528 glomeruli, GSE30529 tubules) + Fig 8 panels
.libPaths(c("~/software/R_libs", .libPaths()))
suppressPackageStartupMessages({
  library(rms); library(pROC); library(ggplot2); library(reshape2)
})
base <- "PROJECT_ROOT"
outdir <- file.path(base, "output", "step8")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

mg <- readLines(file.path(base, "output", "step6", "model_genes.txt"))
coefs <- read.csv(file.path(base, "output", "step6", "LASSO_model_genes.csv"), stringsAsFactors = FALSE)
cf <- coefs$coefficient; names(cf) <- coefs$gene

load_ds <- function(gse) {
  expr <- read.csv(file.path(base, "data", paste0(gse, "_norm_expr.csv")),
                   check.names = FALSE, stringsAsFactors = FALSE)
  rownames(expr) <- expr$gene; expr$gene <- NULL; expr <- as.matrix(expr)
  pheno <- read.csv(file.path(base, "data", paste0(gse, "_pheno.csv")), stringsAsFactors = FALSE)
  list(expr = expr, pheno = pheno)
}

dca_nb <- function(pred, yy, thr = seq(0.01, 0.99, 0.01)) {
  n <- length(yy); prev <- mean(yy)
  nb_m <- sapply(thr, function(t) {
    tp <- sum(pred >= t & yy == 1); fp <- sum(pred >= t & yy == 0)
    (tp - fp * (t / (1 - t))) / n
  })
  nb_all <- prev - (1 - prev) * (thr / (1 - thr))
  list(thr = thr, model = nb_m, all = nb_all, none = rep(0, length(thr)))
}

auc_rows <- data.frame()

# 6-gene training-cohort model (MIR644A absent on GPL571) used for validation nomograms;
# its coefficients are stable, unlike within-cohort refits (LPL shrinks to ~0 in GSE30529)
train6 <- load_ds("GSE96804")
ytr6 <- factor(train6$pheno$group, levels = c("Control", "DN"))
present_tr <- setdiff(mg, "MIR644A")
df_tr6 <- as.data.frame(t(train6$expr[present_tr, ])); df_tr6$y <- ytr6
dd6 <- datadist(df_tr6); options(datadist = "dd6")
f6 <- lrm(y ~ ., data = df_tr6)
print(round(f6$coefficients, 4))

for (gse in c("GSE30528", "GSE30529")) {
  cat("\n====== external validation:", gse, "======\n")
  ds <- load_ds(gse)
  expr <- ds$expr; pheno <- ds$pheno
  stopifnot(identical(pheno$gsm, colnames(expr)))
  y <- factor(pheno$group, levels = c("Control", "DKD"))
  present <- mg[mg %in% rownames(expr)]
  missing <- setdiff(mg, present)
  cat("model genes present:", length(present), " missing:", paste(missing, collapse = ","), "\n")
  Xv <- t(expr[present, , drop = FALSE])
  cfv <- cf[present]
  rs <- as.numeric(Xv %*% cfv)
  if (length(missing) > 0) cat("NOTE: RiskScore computed on", length(present),
                               "genes only (missing probes on GPL571)\n")
  write.csv(data.frame(gsm = pheno$gsm, group = pheno$group, RiskScore = rs),
            file.path(outdir, paste0(gse, "_riskscore.csv")), row.names = FALSE)

  # ROC / AUC
  roc_obj <- roc(response = y, predictor = rs, quiet = TRUE)
  ci <- ci.auc(roc_obj)
  cat(sprintf("%s AUC = %.4f (95%%CI %.4f-%.4f)\n", gse, auc(roc_obj), ci[1], ci[3]))
  auc_rows <- rbind(auc_rows, data.frame(dataset = gse, predictor = "RiskScore",
                                         AUC = auc(roc_obj), CI_low = ci[1], CI_high = ci[3]))
  svg(file.path(outdir, paste0("fig7b_roc_", gse, ".svg")), width = 6.5, height = 6)
  plot(roc_obj, col = "#B2182B", lwd = 2.2,
       main = paste0("ROC in ", gse, ifelse(gse == "GSE30528", " (glomeruli)", " (tubules)")),
       print.auc = TRUE, cex.main = 0.95, cex.lab = 0.95)
  grid()
  dev.off()

  # nomogram: 6-gene training-cohort model (stable coefficients) applied to this cohort
  options(datadist = "dd6")
  svg(file.path(outdir, paste0("fig7a_nomogram_", gse, ".svg")), width = 12, height = 7)
  plot(nomogram(f6, fun = plogis, funlabel = "P(DKD)",
       fun.at = c(0.01, 0.5, 0.99), nint = 4),
       main = paste0("Nomogram (6-gene training-cohort model): ", gse), cex.axis = 0.7)
  dev.off()
  # within-cohort refit only for calibration/DCA
  df <- as.data.frame(Xv[, present]); df$y <- y
  dd <- datadist(df); options(datadist = "dd")
  fv <- lrm(y ~ ., data = df, x = TRUE, y = TRUE)
  pred <- predict(fv, type = "fitted")

  # calibration (apparent for small n; boot fallback)
  cal_ok <- try({
    cal <- calibrate(fv, method = "boot", B = 200)
    svg(file.path(outdir, paste0("fig7c_calibration_", gse, ".svg")), width = 7, height = 7)
    par(mar = c(9.5, 4.5, 3, 2))
    plot(cal, xlab = "Predicted probability of DKD", ylab = "Observed proportion",
         main = paste0("Calibration (bootstrap) in ", gse), cex.main = 0.95)
    abline(0, 1, lty = 2, col = "grey50")
    dev.off()
    TRUE
  }, silent = TRUE)
  if (!inherits(cal_ok, "try-error")) {
    cat("calibration boot ok\n")
  } else {
    cat("calibration boot failed, apparent:\n")
    svg(file.path(outdir, paste0("fig7c_calibration_", gse, ".svg")), width = 7, height = 7)
    cal <- calibrate(fv, method = "apparent")
    par(mar = c(9.5, 4.5, 3, 2))
    plot(cal, xlab = "Predicted probability of DKD", ylab = "Observed proportion",
         main = paste0("Calibration (apparent) in ", gse), cex.main = 0.95)
    abline(0, 1, lty = 2, col = "grey50")
    dev.off()
  }

  # DCA
  dca <- dca_nb(pred, as.integer(y == "DKD"))
  ddf <- data.frame(threshold = rep(dca$thr, 3),
                    NB = c(dca$model, dca$all, dca$none),
                    Strategy = rep(c("Model", "Treat all", "Treat none"), each = length(dca$thr)))
  write.csv(data.frame(threshold = dca$thr, NB_model = dca$model, NB_all = dca$all, NB_none = 0),
            file.path(outdir, paste0(gse, "_DCA_netbenefit.csv")), row.names = FALSE)
  pD <- ggplot(ddf, aes(x = threshold, y = NB, color = Strategy)) +
    geom_line(linewidth = 0.9) + coord_cartesian(ylim = c(-0.05, max(0.3, max(dca$model) * 1.05))) +
    labs(title = paste0("Decision curve analysis in ", gse),
         x = "Threshold probability", y = "Net benefit") +
    theme_bw(base_size = 10)
  ggsave(file.path(outdir, paste0("fig7d_dca_", gse, ".svg")), pD, width = 7, height = 5,
         device = grDevices::svg)

  # single-gene ROC AUCs (Fig 8B equivalent)
  for (g in present) {
    r <- roc(response = y, predictor = Xv[, g], quiet = TRUE)
    c2 <- ci.auc(r)
    auc_rows <- rbind(auc_rows, data.frame(dataset = gse, predictor = g,
                                           AUC = auc(r), CI_low = c2[1], CI_high = c2[3]))
  }
  # model gene expression boxplots
  bx <- melt(expr[present, , drop = FALSE])
  colnames(bx) <- c("gene", "gsm", "expr")
  bx$group <- y[match(bx$gsm, pheno$gsm)]
  pw <- sapply(present, function(g) wilcox.test(expr[g, ] ~ y)$p.value)
  lab <- data.frame(gene = present, p = signif(pw, 3), x = 1.4, y = max(bx$expr) * 0.95)
  p8 <- ggplot(bx, aes(x = group, y = expr, fill = group)) +
    geom_boxplot(outlier.size = 0.7) + facet_wrap(~gene, scales = "free_y") +
    geom_text(data = lab, aes(x = x, y = y, label = paste0("p=", p)),
              inherit.aes = FALSE, size = 2.6) +
    scale_fill_manual(values = c(Control = "#56B4E9", DKD = "#E69F00"), guide = "none") +
    labs(title = paste0("Model genes: DKD vs Control in ", gse), x = NULL, y = "Expression") +
    theme_bw(base_size = 9) + theme(strip.text = element_text(face = "italic"))
  ggsave(file.path(outdir, paste0("fig8a_modelgene_box_", gse, ".svg")), p8,
         width = 10, height = 6, device = grDevices::svg)
}

# ---- Fig 8 (training set): single-gene ROC curves + correlation matrix ----
train <- load_ds("GSE96804")
ytr <- factor(train$pheno$group, levels = c("Control", "DN"))
Xtr <- t(train$expr[mg, ])
svg(file.path(outdir, "fig8b_single_gene_roc_training.svg"), width = 6.5, height = 6)
cols <- rainbow(length(mg)); ltyv <- rep(1:2, length.out = length(mg))
plot(x = NULL, y = NULL, xlim = c(1, 0), ylim = c(0, 1),
     xlab = "1 - Specificity", ylab = "Sensitivity",
     main = "Single-gene ROC curves (training set GSE96804)", cex.main = 0.95)
leg_lab <- character(0); leg_col <- character(0); leg_lty <- numeric(0)
for (i in seq_along(mg)) {
  r <- roc(response = ytr, predictor = Xtr[, mg[i]], quiet = TRUE)
  lines(r, col = cols[i], lwd = 1.8, lty = ltyv[i])
  ci2 <- ci.auc(r)
  leg_lab <- c(leg_lab, sprintf("%s (AUC=%.3f)", mg[i], auc(r)))
  leg_col <- c(leg_col, cols[i]); leg_lty <- c(leg_lty, ltyv[i])
}
legend("bottomright", legend = leg_lab, col = leg_col, lwd = 1.8, lty = leg_lty,
       cex = 0.72, bty = "n")
dev.off()
for (g in mg) {
  r <- roc(response = ytr, predictor = Xtr[, g], quiet = TRUE)
  c2 <- ci.auc(r)
  auc_rows <- rbind(auc_rows, data.frame(dataset = "GSE96804_training", predictor = g,
                                         AUC = auc(r), CI_low = c2[1], CI_high = c2[3]))
}
write.csv(auc_rows, file.path(outdir, "AUC_all_predictors.csv"), row.names = FALSE)
cat("\nAUC table:\n"); print(auc_rows, digits = 3)

# correlation matrix of model genes (training set)
cm <- cor(Xtr)
svg(file.path(outdir, "fig8c_correlation_matrix.svg"), width = 7, height = 6)
par(mar = c(6, 6, 3, 1))
image(x = seq_len(ncol(cm)), y = seq_len(nrow(cm)), z = t(cm)[, nrow(cm):1],
      col = colorRampPalette(c("#2166AC", "white", "#B2182B"))(100),
      xaxt = "n", yaxt = "n", xlab = NULL, ylab = NULL,
      main = "Model gene correlation (training set)")
axis(1, at = seq_len(ncol(cm)), labels = mg, las = 2, cex.axis = 0.8)
axis(2, at = seq_len(nrow(cm)), labels = rev(mg), las = 2, cex.axis = 0.8)
for (i in seq_len(ncol(cm))) for (j in seq_len(nrow(cm))) {
  v <- cm[j, i]
  text(i, nrow(cm) - j + 1, sprintf("%.2f", v), cex = 0.7,
       col = ifelse(abs(v) > 0.5, "white", "black"))
}
box()
dev.off()
cmel <- melt(cm); colnames(cmel) <- c("g1", "g2", "r")
write.csv(cmel, file.path(outdir, "model_gene_correlation.csv"), row.names = FALSE)
cat("\nSTEP8 DONE\n")


############################################################
# FILE: step9_immune.py
############################################################
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Step 9: Immune infiltration on GSE96804 with LM22-derived signatures (ssGSEA).
Paper used CIBERSORT+LM22 (SVR); CIBERSORTx requires a token, so we use the
same LM22 gene sets (547 marker genes x 22 cell types, mastR/CIBERSORT source)
with ssGSEA deconvolution. Outputs Fig 9A/B/C equivalents.
"""
import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.stats import spearmanr
import gseapy

BASE = "PROJECT_ROOT"
OUT = os.path.join(BASE, "output", "step9")
os.makedirs(OUT, exist_ok=True)

expr = pd.read_csv(os.path.join(BASE, "data", "GSE96804_norm_expr.csv"), index_col=0)
pheno = pd.read_csv(os.path.join(BASE, "data", "GSE96804_pheno.csv"))
expr.columns = [c.strip() for c in expr.columns]
assert list(expr.columns) == list(pheno["gsm"])

lm22 = pd.read_csv(os.path.join(BASE, "data", "LM22_signature.csv"), index_col=0)
cell_types = list(lm22.columns)
print("LM22 signature:", lm22.shape, " cell types:", len(cell_types))

gene_sets = {ct: sorted(set(lm22.index[lm22[ct] == 1])) for ct in cell_types}
for ct, g in gene_sets.items():
    print(ct, len(g))

# ssGSEA on log2-scale expression (genes x samples)
res = gseapy.ssgsea(data=expr, gene_sets=gene_sets, outdir=None,
                    sample_norm_method="rank", min_size=5, no_plot=True)
r2 = res.res2d
print("res2d columns:", list(r2.columns))
if {"Name", "Term", "NES"}.issubset(r2.columns):
    nes = r2.pivot(index="Name", columns="Term", values="NES")
else:  # older layout: rows = gene sets, cols = samples
    nes = r2.T
nes.index = [str(i).strip() for i in nes.index]
nes.columns = [str(c).strip() for c in nes.columns]
# align to sample order
nes = nes.loc[[c for c in expr.columns if c in nes.index]]
print("ssGSEA NES matrix:", nes.shape)
nes.to_csv(os.path.join(OUT, "immune_nes_matrix.csv"))

# ---- Fig 9A: stacked bar of per-sample immune composition ----
pos = nes.clip(lower=0)
frac = pos.div(pos.sum(axis=1), axis=0)
fig, ax = plt.subplots(figsize=(14, 5.5))
bottom = np.zeros(len(frac))
cmap = plt.get_cmap("tab20")
order = np.argsort(pheno["group"].values, kind="stable")
for i, ct in enumerate(cell_types):
    v = frac[ct].values[order]
    ax.bar(range(len(frac)), v, bottom=bottom, color=cmap(i % 20), width=0.8,
           label=ct, edgecolor="none")
    bottom = bottom + v
ax.set_xticks(range(len(frac)))
ax.set_xticklabels([])
ax.set_xlabel("Samples (Control on the left, DN on the right; n = 61)", fontsize=9)
ax.set_ylabel("Relative abundance\n(ssGSEA NES, positive-normalized)", fontsize=9)
ax.set_title("Immune cell composition per sample (GSE96804, LM22 signatures)")
ax.legend(bbox_to_anchor=(1.01, 1), loc="upper left", fontsize=7.5, ncol=1)
plt.tight_layout()
plt.savefig(os.path.join(OUT, "fig9a_immune_stacked.svg"), bbox_inches="tight")
plt.close()

# ---- Fig 9B: Spearman correlation heatmap among 22 cell types ----
rho = nes.corr(method="spearman")
rho.to_csv(os.path.join(OUT, "immune_celltype_spearman.csv"))
fig, ax = plt.subplots(figsize=(9.5, 8))
im = ax.imshow(rho.values, cmap="RdBu_r", vmin=-1, vmax=1)
ax.set_xticks(range(len(cell_types)))
ax.set_yticks(range(len(cell_types)))
ax.set_xticklabels(cell_types, rotation=90, fontsize=6)
ax.set_yticklabels(cell_types, fontsize=6)
for i in range(len(cell_types)):
    for j in range(len(cell_types)):
        v = rho.values[i, j]
        if i != j:
            ax.text(j, i, f"{v:.2f}", ha="center", va="center",
                    fontsize=4.5, color="white" if abs(v) > 0.6 else "black")
fig.colorbar(im, shrink=0.8)
ax.set_title("Spearman correlation among 22 immune cell types (ssGSEA NES)")
plt.tight_layout()
plt.savefig(os.path.join(OUT, "fig9b_celltype_corr_heatmap.svg"), bbox_inches="tight")
plt.close()
# strongest positive/negative pairs
pairs = []
for i in range(len(cell_types)):
    for j in range(i + 1, len(cell_types)):
        pairs.append((cell_types[i], cell_types[j], rho.values[i, j]))
pairs = sorted(pairs, key=lambda x: x[2])
pd.DataFrame(pairs, columns=["cellA", "cellB", "spearman_r"]).to_csv(
    os.path.join(OUT, "immune_celltype_pairs.csv"), index=False)
print("strongest negative:", pairs[:2])
print("strongest positive:", pairs[-2:])

# ---- Fig 9C: model gene vs immune cell Spearman correlation bubble plot ----
mg = [g.strip() for g in open(os.path.join(BASE, "output", "step6", "model_genes.txt"))]
mg = [g for g in mg if g in expr.index]
corr = pd.DataFrame(index=mg, columns=cell_types, dtype=float)
pval = pd.DataFrame(index=mg, columns=cell_types, dtype=float)
for g in mg:
    for ct in cell_types:
        r, p = spearmanr(expr.loc[g, nes.index].values, nes[ct].values)
        corr.loc[g, ct] = r
        pval.loc[g, ct] = p
corr.to_csv(os.path.join(OUT, "model_gene_immune_corr.csv"))
pval.to_csv(os.path.join(OUT, "model_gene_immune_pval.csv"))

sig = (pval < 0.05).values
fig, ax = plt.subplots(figsize=(11, 4.2))
for i, g in enumerate(mg):
    for j, ct in enumerate(cell_types):
        r = corr.values[i, j]
        ax.scatter(j, i, s=60 * abs(r) + 4,
                   c="#B2182B" if r > 0 else "#2166AC",
                   alpha=0.9 if sig[i, j] else 0.25,
                   edgecolors="none")
ax.set_yticks(range(len(mg)))
ax.set_yticklabels(mg, fontsize=8, style="italic")
ax.set_xticks(range(len(cell_types)))
ax.set_xticklabels(cell_types, rotation=90, fontsize=6)
ax.invert_yaxis()
ax.set_title("Spearman correlation: model genes vs immune cells (GSE96804)\n"
             "filled = p<0.05; red = positive, blue = negative; size ~ |r|")
mappable = plt.cm.ScalarMappable(cmap="RdBu_r")
plt.tight_layout()
plt.savefig(os.path.join(OUT, "fig9c_gene_celltype_bubble.svg"), bbox_inches="tight")
plt.close()
print("step9 done")
